WTreucy 
4 ’ Tas 
aA 


+ 


PIATTE RA A CTRY 
pS 


) 4 . ° iv} 
— pe Ge ar PS es 


EDITED BY 


T..& FRY TH. v. KARMAN 
W. PRAGER I. S. SOKOLNIKOFF 
J. L. SYNGE 


ITH THE COLLAVORATION OF 


L. N. BRILLOUIN 
W. FELLER 

J. J. N. GOODIER 

P. LE CORBEILLER F, D. MURNAGHAN 

W. R. SEARS R. V. SOUTHWELL 

3. P. TIMOSHENKO H. S. TSIEN 





JANUARY + 1946 NuMBER 4 


























QUARTERLY 


OF 


APPLIED MATHEMATICS 


H. L. DRYDEN 
J. M. LESSELLS 


H. BATEMAN 

J. P. DEN HARTOG 

K. O. FRIEDRICHS 

G. E. HAY 

S. A. SCHELKUNOFF 
SIR GEOFFREY TAYLOR 


EDITED BY 


1. Re 
W. PRAGER 
J. L. SYNGE 


WITH THE COLLABORATION OF 


M. A. BIOT 

H. W. EMMONS 

J. A. GOFF 

P. LE CORBEILLER 
W. R. SEARS 

S. P. TIMOSHENKO 


TH. v. KARMAN 
I. S. SOKOLNIKOFF 


L. N. BRILLOUIN 
W. FELLER 

J. N. GOODIER 

F. D. MURNAGHAN 
R. V. SOUTHWELL 
Hi. S. TSIEN 





VotumeE III 











(,EORGI 


Printed by the 


BANTA 


Menasha 


»~UBLISHING COMPANY 


Wisconsin 


CONTENTS 


T. Alfrey: Methods of representing the properties of viscoelastic materials . 

H. J. Barten: On the deflection of a cantilever beam ee eo” 

K. E. Bisshopp: The inverse of a stiffness matrix... awe 

K. E. Bisshopp and D. C. Drucker: Large deflection of cantileves benme 5» We 

D. R. Blaskett and H. Schwerdtfeger: A formula for the solution of an arbitrary 
analytic equation ‘ wineries teria tur? oll sill. ofnegeag 

M. H. Blewett: (See //. Poritsky) 

E. Bromberg and J. J. Stoker: Non-linear theory of curved elastic sheets . 

G. F. Carrier: On the non-linear vibration problem of the elastic string . 

G. F. Carrier: On the vibrations of the rotating ring. . . é 

P. Y. Chou: On velocity correlations and the solutions of the equations of 

turbulent fluctuation ; er 

P. Y. Chou: Pressure flow of a turbulent fluid between two infinite pataliel 

planes Low eS & oo oes Se, Le Se eee eee 
N. Coburn: The Karman-Tsien pressure-volume relation in the two-dimen- 
sional supersonic flow of compressible fluids 


J. W. Craggs and C. J. Tranter: The capacity of twin cable... . 268, 

C. H. Dix, C. Y. Fu and E. W. McLemore: Rayleigh waves and ie surface 
reflections... -. i f & .4« ce 

D. C. Drucker: (See K. E. " Bisshopp) 

C. M. Fowler: Analysis of numerical solutions of transient heat-flow problems . 

T. C. Fry: Some numerical methods for locating roots of polynomials 

C. ¥. Fu: (See C. H. Déx) 

W. Horenstein: On certain integrals in the theory of heat conductions 


J. C. Jaeger: Diffusion in turbulent fiow between parallel planes 
Th. v. Karman and H. S. Tsien: Lifting-line theory for a wing in non- naiboren 
flow Seek eS Oe a ee eae ee ae 
R. King and D. Middleton: The cylindrical antenna; current and impedance . 
E. G. Kogbetliantz: Quantitative interpretation of maps of magnetic and 
gravitational anomalies by mathematical methods 
W. Kohn: The spherical gyrocompass 
C. C. Lin: On the stability of two- Menenabeinil parallel Sows. 
I. General theory 
Il. Stability in an inviscid fluid 
[1]. Stability in a viscous fluid a eee ee 
A. N. Lowan: On the problem of heat conduction in a semi- inGnite radiating 
wire ° . ° ° 
E. W. McLemore: (See C. HT. Dix) 
D. Middleton: (See R. King) 
I. Opatowski: Cantilever beams of uniform strength , 
A. A. Popoff: A new method of integration by means of orthagenalits foci , 
H. Poritsky and M. H. Blewett: A method of solution of field problems by 
means of overlapping regions... ae ae i 
W. Prager: On plane clastic strain in doubly commectedl doen Lins 


143 
275 

82 
272 
266 
246 
157 
235 

38 
198 


106 
380 


361 
89 


183 
210 


302 


Coo ul 
“su 


76 
166 


336 
377 











A. Preisman: Graphical analyses of nonlinear circuits . . . 

M. Richardson: The pressure distribution on a body in shear flow 

H. A. Robinson: On A. A. Popoff’s method of integration by orthogonality foci . 

S. A. Schaaf: A cylinder cooling problem ek le Cr eee 

S. A. Schelkunoff: Solution of linear and slightly nonlinear differential equa- 
tions 


H. Schwerdtfeger: (See D. R. Blaskett) 

C. H. W. Sedgewick: On plastic bodies with rotational symmetry 

J. J. Stoker: (See E. Bromberg) 

S. Timoshenko: On the treatment of discontinuities in beam deflection prob- 


ia ok ee 
C. J. Tranter: (See J. W. Craggs) 
H. S. Tsien: (See Th. v. Kérmén) 
A. Vazsonyi: On rotational gas flows 2 Oe ee ee ee eee 
S. E. Warschawski: On Theodorsen’s method of conformal mapping of nearly 
circular regions 
Book Reviews 


Bibliographical lists . 184, 27 


85 
75 





277 


QUARTERLY OF APPLIED MATHEMATICS 





Vol. III JANUARY, 1946 No. 4 





ON THE STABILITY OF TWO-DIMENSIONAL 
PARALLEL FLOWS 


PART III.—STABILITY IN A VISCOUS FLUID* 


BY 


cc. <. Lan 
Guggenheim Laboratory, California Institute of Technology 


11. General considerations. The investigations in Part II* of the stability char- 
acteristics of a parallel flow in an inviscid fluid led to very useful information. In 
the first place, they enable us to visualize the effect of pressure forces very clearly. 
In the second place, the results can be used as a guide for studying the stability prob- 
lem in a real fluid, since instability is expected to occur only for sufficiently large 
Reynolds numbers. Thus, if we know the general characteristics of “inviscid stability” 
for a given velocity distribution, some stability characteristics in a real fluid can be 
obtained by considering a modification of these results by the effect of viscosity. Such 
a consideration was first made by Heisenberg, [14]+ who demonstrated that the ef- 
fect of viscosity is generally destabilizing at very large Reynolds numbers. There are, 
however, a few points to be supplemented in his discussion. We shall therefore study 
this problem in some detail in §12. 

To get a good understanding of the stability problem, we want to know the fol- 
lowing points for any given class of velocity distribution. First of all, we want to 
know whether this class of flows is stable for all Reynolds numbers. Secondly, if it is 
stable for certain Reynolds numbers with respect to disturbances of certain wave- 
lengths, but unstable under other conditions, we want to know the general nature of 
the curve c;(a, R) =0 which separates regions of stability and instability in the a, R 
plane. Thirdly, such a curve will be expected to show a minimum in R, below which 
all small disturbances are damped out. It is therefore desirable to be able to calculate 
this minimum critical Reynolds number rapidly. 

In the next section, we shall solve these problems for two classes of flows; namely, 
(a) velocity distributions of the symmetrical type, and (b) velocity distributions of 
the boundary-layer type. Indeed, it will be shown that in these cases the flow is always 
unstable for sufficiently large Reynolds numbers, whether the velocity curve has a point 
of inflection or not. The curve of neutral stability c;(a, R) =0 will be shown to belong 
to either of the types shown in Fig. 9. When the velocity curve has no point of inflec- 
tion, the two asymptotic branches of the curve have the common asymptote a=0 


* Received July 18, 1945. Parts I and II of this paper appeared in this Quarterly 3, 117-142, and 
218-234 (1945). 

** Now at Brown University. 

Tt The figures in brackets refer to titles in the Bibliography at the end of Part I. 
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(Fig. 9a). When there is a point of inflection, one branch has the asymptote a=0, 
while the other has the asymptote a=a,>0 (Fig. 9b). In either case, it will be shown 
that the region inside the loop is the region of instability. It should be observed how 
these results fit in with those of Rayleigh and Tollmien for the inviscid case, and 
with the destabilizing effect of viscosity noted by Heisenberg for large Reynolds 
numbers. Simple formulae will be derived to express the asymptotic branches of the 
curves in Fig. 9. In fact, it is by means of these asymptotic behaviours and a criterion 
of stability of Synge [63], that the results mentioned above are established. Simple 
rules will also be given, by which the minimum critical Reynolds number can be very 
easily obtained from quantities involving very simple integral and differential ex- 
pressions of the velocity distribution w(y). Very little numerical labor is involved for 
the calculation in any particular case. 

Heisenberg also discussed the general shape of the curve c;(a, R) =0. However, 
his argument does not appear to be very decisive, and some of his results are not well 
stated. He did not obtain the asymptotic forms of the a(R) curve as given below. He 
also tried to estimate the order of magnitude of the critical Reynolds number, but 
did not try to make an approximate calculation in terms of simple differential and 
integral expressions of w(y) (loc. cit., p. 600). 

In order to obtain definite numerical results, which may be subjected to experi- 
mental verification, we shall apply our theory to the stability problem of special 
velocity distributions. In §13, we shall give the results of calculation of the neutral 
curves in (a) the Blasius case and (b) the Poiseuille case. The method of calculation 
and its numerical accuracy will be given in the Appendix. Frequent reference to the 
equations in it will therefore be made in the following sections. 

12. Heisenberg’s criterion and the general characteristics of the curve of neu- 
tral stability. We shall now proceed to study the general stability characteristics in 
a viscous fluid as indicated above. In the first place, we shall develop Heisenberg’s 
criterion in a slightly improved form. We shall then restrict ourselves to velocity dis- 
tributions of the symmetrical type and of the boundary-layer type. For these cases, 
we shall prove the results summarized in the next section 

a) Heisenberg’s criterion. Along the neutral curve 


ci(a, R) = 0 


(if it exists), all the parameters a, c, and R are functions of one of them, say, R. Let us 
restrict ourselves to cases where a and ¢ do not approach zero along the neutral curve 
as R becomes infinite. Then, the approximations (6.27) are valid for sufficiently large 
values of R. By using (6.26), (6.24) and (6.27), we can transform (6.13) into the fol- 
lowing: 


Pa e~ri/4 20 20 
0? "Kons (c) = —— ; a?" Kon(c) + we (1 = Cc) a?"Ko, @} 
~ ne) © ek — 2 ~ ™ 
etil4 s 20 00 
. bs a?" H5,(c) + wc > a" Kmnin(Ct ° 
Vv (aRc)5 \ n=0 n=( 


If there is an inviscid neutral disturbance with c=c,, a=a,, we have 


~ 9 
> as Kensi(cs) = 0. 


n=0 
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Hence, when a=a,, and c=c,+Ac, differing but slightly from c,, we have 








L 2 on Kon(Cs) - ¥: a, Hon (Cs) 
ee orn. oe n=0 Ore n=0 
" Vv a,R(1 — ¢,)5 ~ 2n bi Va,Re - 2n , , si 
z as K on+1(Cs) ) as Kons 1(Cs) 
n=0 n=0 


Similar considerations of (6.14), (6.15), and (6.17) give respectively 


oo 














2n 
isa Qs Hon+1(Cs) 
eri pi 
is 2 -S . ’ (12.2) 
Vv a,Re5 os 2n al 
) ® As Kn+2(Cs) 
n=0 
a 2n 
- De & Hon(cs) 
n= , (12.3) 
V as C, = 2n alt 
z. As Kon+1(Cs) 
n=0 
a 2n 2 = 2n 
Wa bs as Aon-1 (s) + (1 — Cs) >» as "Hon( Cs) 
a = a ad 
dc = Ta = _ : (12.4) 
lt, Qn 2 2n+1 7 
7. as Kan(Ce) + (1 ts Cs) ps as ‘ Kon+1(Ce) 
n=1 n=0 


In general, it is not very easy to determine whether Ac;>0 or <0. However, when 
c, and @, are both small, but not zero, all the above equations will reduce to 


Ac = e*/4/\/a,Re® Ki (c.), (12.5) 
after we make use of the reductions corresponding to (1) and (2) of the Appendix. 
Now, when ¢ is small, 

1 we’ 


wic wy, 3 





(log ¢ + ir) + O(1); 





v2 
K,(c) -f dy(w—c)? = — 
yy; 

hence, it can be easily verified that KY (c) is approximately real and positive for small 
real values of c. Hence, (12.5) shows that Ac;>0 in every case. The disturbance with 
wave length 27/a,, neutrally stable in the inviscid case, is unstable when viscous 
forces are considered. This result was first obtained by Heisenberg, and may be 
formulated as follows: 

HEISENBERG’S CRITERION. If a velocity profile has an “inviscid” neutral disturbance 
with non-vanishing wave number and phase velocity, the disturbance with the same wave 
number is unstable in the real fluid when the Reynolds number is sufficiently large. 

In Heisenberg’s original discussion, only the first type of motion is considered. 
The last equation on p. 597 of his paper is essentially our equation (12.1) with all 
the terms in a? dropped. Evidently, the above arguments hold only for a,, c, #0. It will 
be seen later from Fig. 9 that one neutral disturbance with a=0 and c=0 for infinite 
R, is actually stable for finite values of R. 
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b) Asymptotic behavior of the neutral a(R) curve. We shall now study the general 
asymptotic behavior of the neutral curve, assuming that it exists. The answer to the 
existence problem will be evident during the course of the investigation. For large 
values of aR, we would generally expect z of (6.28) to be much greater than 1, but it 
is also possible for z to approach a finite value or zero. We shall therefore discuss both 
possibilities. 

For large values of 2, we have approximately, 


7, =1, F, = 1/22? = wl /V/2aRe3, (12.6) 
where 7; is small. If we refer to (5) and (7) of Appendix, we see that the imaginary 
part v of the right-hand side member of those equations depends on that of wy cK,(c) 
and those of the integrals Hi, Hz, M’s and N’s. If a and ¢ are small, which will be 
verified a posteriori, we have only the contribution from the first term; thus,* 

ww’ 


v= — TW — 
w 





for w=c. (12.7) 
‘3 
By using (12.6) and (12.7), we see that the equations (8) in the Appendix can be 


approximated by 
“w= 1, (12.8) 


and 
,I0=-—- -— CWo = wy / Vv 2aRe?, (12.9) 
Wi 
These are the equations for determining a relation a(R), if we eliminate c between 
them. 

In the case where aRc approaches a finite limit as aR— ©, c must approach zero. 
Hence, v must approach zero, and from (8) of Appendix, 7; must also approach zero. 
From the curve for 7;(z), we see that this happens for z = 2.294, for which 7, = 2.292. 
Then, using (9) of Appendix, (12.6), and (12.7), we have 


c3, g = 2.294, (12.10) 


and 
u = F, = 2.292. (12.11) 


The two types of relations (12.8), (12.9) and (12.10), (12.11) evidently correspond 
to two different branches of the a(R) curve (cf. Fig. 9). These conditions can be 
satisfied in cases (2a) and (3) of section 6 (cf. (11.5), (11.7) of Appendix), but it 
appears to be difficult in case (2b) (cf. (11.6) of Appendix). 

CAsE (2a). Symmetrical velocity distribution with symmetrical $(y). We consider 
the cases where both a and ¢ are small. By using (12.5) and noting that u takes on a 
finite value in either (12.8) or (12.11), we see that we must have 

wic A) 
u = ——», where Hy = A,(0) = f wdy; (12.12) 
Hypa? Vi; 

* In fact, the other terms never give considerable contributions to the imaginary part even for only 

moderately small values of a and c. This point will be discussed in the Appendix. The approximation (12.7) 





will be used for all later calculations. 
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i.e., C must approach zero as fast as a?. The asymptotic behavior of the a(R) curves 
as given by (12.8)—(12.11) are as follows: 


R = (wi /22*Hiows da", ¢ = (Hy/w;)e, (first branch), (12.13) 


R 


wi'(2/F,Hi») a’, c = (Hw7,/w)e2, (second branch), (12.14) 


where ¥,=2.292, z=2.294 approximately. 
Case (3). Boundary-layer profile. Here, the equation corresponding to (12.12) is 


(cf. (7) of Appendix) 
u= w,c/a, (12.15) 


i.e., © must approach zero as fast as a. Note that in the previous case, the relation 
(12.12) depends both on wy and on /¥w*dy. Here, it depends only on the initial slope 
of the velocity curve w/ . The two branches of the a(R) curve for large values of R are 


R = (wi "/2n?wo '*) a, c=a/w, (first branch), (12.16) 


and 


I 


R = w (2/F,)*a, c= aF,/wi (second branch), (12.17) 


where 7, = 2.292, z=2.294 approximately. 

Effect of varying curvature in the curve of velocity distribution. In either case, the 
second branch of our asymptotic curve depends very little upon the shape of the 
velocity profile, while the first branch depends very much upon it through the term 
wg’. This fact will enable us to correlate our present results with the inviscid inves- 
tigations of Rayleigh and Tollmien, as discussed in Part II. 

In all the cases considered, we have w’’ <0 for y<ye2 but sufficiently near to it. 
If w’’(y) never vanishes for y;<y<ye, we can replace wg’ by wi’ in the expressions 
(12.13) and (12.16). In general, 





Now, for a flow which is essentially parallel, the boundary condition AAy =0, which 
holds on the solid wall for all two-dimensional laminar flows, can be easily verified 
to be equivalent to w/" =0. Hence, we have 


wiv 


YS a 





Ul 
W =w + 
/ 
2w;? 


Thus, if w{/’ =0, but w’’ does not vanish for y1$<y<ye, we have 

R = {2w{" ‘n° H},(w'’)* a, for case (2a), (12.18) 
and 

R = {2w{*/x*(w)*}a-™, for case (3). (12.19) 
In case the velocity profile shows a point of inflection, 


wo =0 for c=. 
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Then, we have approximately 


wo = (wi'/2wi*)(c? — c?). (12.20) 


It can now be seen from (12.13) and (12.16) that R becomes infinite as c approaches 
c,. Let the corresponding value of a be denoted by a,. Then instead of (12.18) and 
(12.19), the following relations hold: 


19 2779 2 —11 4 4, —2 2 
R = }2w{/2°Hy(w;) fa (a* — at)’, at = wic,/Hy, for case (2a), (12.21) 
19 2 r\ 2 —6 2 2, —2 / 
R = {2wi"’/x"(w,")*}a “(a — ai), a, = wic,, for case (3). (12.22) 
Thus, for either a symmetrical or a boundary-layer distribution with a flex, we have 
R~ (a — a,)~? (12.23) 


as R>~, a—a,, c—c,. In all these approximations, we assume a, and c, to be so small 
that the previous approximations still hold, but the qualitative nature of our conclu- 
sions cannot be changed for moderate values of a, and ¢,. 

The general characteristics obtained from our foregoing discussions are summa- 
rized in Table II, and are indicated by the asymptotic branches of the solid lines in 
Fig. 9. Let us proceed to discuss their significance. 

i) It may be expected that the region between the two asymptotic branches of the 
curves represents a region of instability. Thus, every symmetrical or boundary-layer 
profile is unstable for sufficiently large values of the Reynolds number. This point will be 
substantiated below. 

ii) In the cases where w{’ >0, the two branches of curves approach the two differ- 
ent asymptotes a=0 and a=a,, leaving a finite instable region for infinite Reynolds 
number. In the other cases, the two branches approach the same asymptote a=0, 
leaving only the possibility of a neutral disturbance of infinite wave-length at infinite 
Reynolds number. These results agree with Heisenberg’s criterion and the results ob- 
tained from the considerations of an inviscid fluid in Part II of this work. It thus ap- 
pears that the inviscid disturbance with a=0, c=0 is actually not as trivial as it 
may first appear to be, for it is actually the limiting case of neutral disturbances in a 


real fluid. 

c) Existence of self-excited disturbances.* To establish the actual existence of self- 
excited disturbances, we try to show that c;>0 at least in the neighborhood of the 
neutral curve. Indeed, we may regard ¢ as a function of the two independent param- 
eters a and R’ =aR, and show that (0c;/0R’), <0 for the first branch of the curve. This 
is analogous to Heisenberg’s criterion, and demonstrates the same general conclusion 
that the effect of viscosity is destabilizing at large Reynolds numbers. To fix our ideas, 
we shall carry out the analysis for the case of symmetrical profiles. The other case 
can be carried out in a similar manner. 

We begin with the equation 

fola,c) dai 

f(a, ©) $4; 
where f2(a, c) and f(a, c) are given by (6.26) and ¢3:/¢3; is given by (6.28). By using 
those relations, we can transform (6.14) into 


* This section was inserted late in 1944 after discussions with Prof. C. L. Pekeris. He mentioned the 
possibility that the neutral curve might be a curve of minimum damping with stable regions on both sides 


of it. See also Schlichting’s calculations [52]. 


(6.14) 
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or, by further using (6.30) and (6.24), 


H(z) —1 = wic >> a*Kal(c) / Do a" Hnr(c). 


n=l n=1 


We now regard was fixed and consider the variation of c with R’ or z, which is a known 
function of c and R’. Taking logarithmic derivatives on both sides, we have 


p a2” K5n(c) > a2" Hon-3(C) 





F(z 1 a= n= c 
artes = —. 
aie da" Koa(c) a?" an—-r(c) 
n=1 n=l 


So far, a, c, and z are arbitrary. On the neutral curve, c and z are real, and we may 
use the relation (12.13) if we are on the first branch of the neutral curve with large 
values of R’. Indeed, for large values of 2, (12.6) gives 


7m 8s 
F(z) — 1 22 
and 
3dz 1 dR’ 3 dc 
Ss “Ss 2 3s 


By using these relations, it can be easily verified that (12.24) leads to 


Dia™*Kin(c) Do a?*Han-1(0) 








dR’ 3dc 1 ond n=l 

> eae * 

2 2¢ c = nas 
> a?" Ko,(c) ym a?" H o,—1(€) 
n=] n=1 


Remembering that c=O(a?) and noting that 


7 
Wo 


1 
K,(c) = -——- (log ¢ + ir) + O(1), 


WC wo? 





where O(1) is real, we can reduce the above equation to 


Ul 


ee alee so} "$form (log ¢ + in) 
2k’ & ( wi? dco wo? 
or , ” 
1 at +<¢ “(ee (log ¢ + in) |. (12. 25) 
2 2¢ dc we * 


For small values of c, the expression in the brackets has a positive real part and a 
negative imaginary part. Hence, (dc/0R’), has a negative imaginary part. This com- 
pletes the proof. Indeed, if w/’ does not vanish, we have 
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(=) Cc 
= — —w R65, 
OR'/ . 3R’ 
(12.26) 


OC; 2rwi' c? 
as: ns —~ R18, 
OR'/ Ow, R’ 


In the above derivation, it should be noted that all approximations are made by 
neglecting small terms of higher orders in comparison with some terms which have 
been retained. Thus, the conclusion of stability or instability would not be altered by 
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those terms of higher orders. 
d) The minimum critical Reynolds number and the minimum critical wave-length. 


Having demonstrated the instability of the symmetrical and the boundary-layer pro- 
files, we want to answer the following questions. First, does there exist a minimum 
critical Reynolds number, below which the flow is stable for disturbances of all wave- 
lengths? If so, can we get an approximate estimate of its magnitude? Secondly, does 
there exist a minimum wave-length of the disturbance (maximum a@) below which the 
flow is stable at all Reynolds numbers? If so, can we get an approximate estimate of its 
magnitude? We shall see that in trying to answer these questions, we can also roughly 
depict the complete a(R) curve, which separates stability from instability. 

The existence of these minimum values can be most conveniently inferred from a 
condition of stability derived by Synge* from energy considerations. His condition 
reads 
(GR)? < (2a? + 1)(4a4+ 1)/a?, gq = max | w’|. (12.27) 
This condition insures the existence of a minimum critical Reynolds number. It per- 
mits a to become infinite only for R-«. But we know from our previous considera- 
tions that a—a, or 0 as R-~. Hence, we would expect that there exists a maximum 
value of a, above which there is always stability. The neutral curve must therefore 
take the general shape shown in Fig. 9. The asymptotic behaviors of the solid curves 
are drawn in qualitative accordance with (12.6), (12.7), and (12.23); the other parts 
of the solid curves are arbitrarily drawn to indicate the general shape of the curve. 
It is evident that the region outside the curve is the region of stability, and the en- 
closed region is the region of instability. Similar conclusions have been reached by 
Heisenbergf but his arguments and results appear to be somewhat obscure. 

Having established the existence of the minimum critical value of R and the 
maximum critical value of a, we proceed to make an estimate of their magnitude. 
We shall see that our theory permits us to give a quite good approximation to the 
minimum value of aR (cf. (12.30)). Since this roughly corresponds to the minimum 
value of R and also to the maximum value of a@ (as will be clear from the individual 
examples given below), we can get a rough estimate of these values by making a rough 
estimate of a corresponding to the minimum value of aR. 

Using the second equation of (8) of Appendix and the approximation (12.7) for v, 
we have approximately 

ww" 
—- (12.28) 


F(z) = v(c) = rwy os 
w’? 





* Synge, [63], eq. (11.23), p. 258. His \ is our a. The condition is originally stated for plane Couette 
and plane Poiseuille motion; but it is easily seen that it holds for a general velocity distribution with 
g=max|w’|. 

+ Heisenberg, loc. cit., p. 601. 
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Fic, 9, General nature of the curve of neutral stability. The dotted curve is curve of 
stability given by Synge. 


If we recall that z is proportional to c(aR)'/*, this equation determines (aR)'/* as a 
function of c. It can then be easily verified that the minimum value of (@R)'/* occurs 


when 
zi (z) = cv'(c). (12.29) 
If the point where this holds is denoted by z = 20, we have approximately from(11.28)* 
Zo 3 
aR = wi?(=), (12.30) 
¢ 


* Cf. Heisenberg, loc. cit., eq. (29a), p. 602. He put so~1. 
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The point 2 is roughly the value where 7;(z) takes its maximum value, because 
(12.29) is approximately ¥/(z)=0, when c is small. The corresponding value of a 
can be approximately obtained by taking 


u = F,(z0), (12.31) 


in accordance with the first equation of (8) of Appendix, where u is given by the real 
part of (5) or (7) of Appendix, as the case may be. 

e) Approximate rules. We now proceed to make some rough approximations in 
order to obtain simple rules, which are convenient for estimating the minimum criti- 
cal Reynolds number. The condition 7/ (zo) =0 gives 29 = 3.21, where 7,(z0) =1.49 and 
F(Z) =0.58. The corresponding value of c can be obtained from the second equation 
of (8) of Appendix. Putting wu =1.5 in accordance with (12.31) and neglecting second 
order terms of A, we have 

(1 — 2A) = Fi(zo) = 0.58, 

i.€., 
— ww; (1 — 2d)(ww’’/w’*) = 0.58, (12. 32) 
where J is defined by (4) of the Appendix. To find the value of c from this equation, 
it is convenient to plot its left-hand side together with w(y) against y, and read the 
value off the latter curve where the former curve gives the value 0.58. The value of 
c so obtained turns out to be very close to its maximum value along the neutral 
curve, and is approximately the value where R is a minimum. 

The values of a and R must be obtained from more rough approximations. With 
a consultation of the values of the integrals H(c), K(c), M(c) and N(c) given in the 
Appendix, we may derive the following reasonable estimates of a: 


y2 

a? = wi c/H, Hi = f w*dy, for symmetrical profiles, (12.33) 
yl 

a= wie, for boundary layer profiles. (12.34) 


These values turn out to be somewhat lower than the accurate values. With an ap- 
proximate allowance for these inaccuracies and taking round numbers, we get the 
following approximate rules for the minimum critical Reynolds number: 








30w!  /Hrots : 
a ———, for symmetrical profiles, (12.35) 
¢ c 
25wy 
R=- oy for boundary layer profiles. (12.36) 
c 


The calculations have been carried out for the Blasius case and the plane Poiseuille 
case. In the first case, the thickness of the boundary layer is taken so that the initial 
slope is w;’=2. It is found that 


R = 5906 for Poiseuille case, \ 


Ré; = 502 ‘for Blasius case. 


The quantity 6; is the displacement thickness 


5; -f (1 — w)dy = 0.28673, 
0G 
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where y is measured from the solid wall. These values for the minimum critical Rey- 
nolds number agree fairly well with those obtained below from more elaborate numer- 
ical calculations. 

When these estimates of the minimum values of R and the corresponding values 
of a (eqs. (12.32)—(12.36)) are combined with the asymptotic behavior of the a(R) 
curves (eqs. (12.16)—(12.17)), the curve of neutral stability in any case can be 
sketched with fair accuracy with very little labor. 

The maximum value of @ on the neutral curve cannot be very well estimated. It 
is usually somewhat higher than the values of a given by (12.33) and (12.34). 

13. Stability characteristics of special velocity distributions. We shall now apply 
our theory to some special cases in order to obtain numerical results comparable 
with experiments. We take (a) the Blasius case as a typical boundary-layer profile, 
and (b) the plane Poiseuille motion as a typical symmetrical profile. In any case, the 
resultant curve of stability limit should have the general shape discussed in the last 
two sections. Only the results will be given here; the method of calculation and its 
accuracy will be discussed in the Appendix. 

a) Stability of plane Poiseuille flow. The velocity distribution of plane Poiseuille 


motion is given by 


w(y) = 2y— y*, with wi =2, H,(0) = 8/15. (13.1) 


TABLE II, Behavior of R(a) Curve for Large Values of R for Velocity Distributions 
with w’’ <0 for the Main Part of the Profile. 








| First branch | 
| Second branch 























w,"’<0 | wi” =0 | w,"">0 
Symmetrical profile a | an | (a—a,) a7 
Boundary-layer profile a6 a= | (a—a,)* | a 
The two branches of the a(R) curve are given by (cf. (12.13), (12.14)) 
R'/3 = 8.44(a?)-"/6, c¢ = 407/15; 
(13.2) 
R13 = §,96(a?)—7/6, c = 0.611a’. 


The numerical results are shown in Table III and Fig. 10. The significance of the 
column s in the table will be explained in the next section. From the figure, we see 
that the minimum critical Reynolds number occurs at R’/* =17.45, or R=5314, agree- 
ing very well with our previous estimation.* 

Earlier results. The stability of plane Poiseuille flow has been attempted by many 
authors. Comparatively recent papers are those of Heisenberg, [14], Noether [36], 
Goldstein [6], Pekeris [39, 40], Synge [64], and Langer [25]. The papers of Goldstein 
and Synge and one of the papers of Pekeris [39] give definite indication of stability 
at sufficiently low Reynolds numbers. Heisenberg’s paper is in general agreement 


* The values given here are somewhat different from those published before [27], because the com- 
putation of Tietjen’s function has been revised. 
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with the present investigations. He gave only a very rough calculation, whose result 


is reproduced in the figure. It seems that his curve is 


R13 = 13.4(a?)-1/6, 


This is different from our present result (13.2) by a numerical factor. It may be noted 


TABLE III. Stability of Plane Poiseuille Flow. 







































































c Zz a | R | Ss | a? 
0 2.294 0 | 00 | 7a | @ - 
0.05 2.363 0.3056 | 13 .64X 105 .6901 | 0.0934 110.91 
0.10 2.448 0.4603 1.243105 | .6544 | 0.2119 49 .90 
0.15 2.540 0.6024 31048 .6192 | 0.3629 31.43 
0.20 2.668 0.7506 12024 .5752 0.5634 22.91 
0.25 2.868 0.9263 6108 .5161 0.8580 18.28 
0.266 3.012 1.0101 5369 -4795 1.0203 | 17.51 
0.270 3.080 1.0414 5314 .4637 1.0845 17.45 
0.272 os | 1.0836 5659 .4358 | 1.1741 17.82 
0.270 3.240 1.0888 } 5920 .4298 1.1854 18.09 
0.266 3.320 1.1007 6602 .4144 | 1.2015 18.76 
0.25 3.495 1.1033 9287 .3836 ‘2073 21.02 
0.20 3.857 1.0254 26597 } .3309 1.0514 29.85 
0.15 4.152 0.8824 92529 .2963 | 0.7787 45 .23 
0.10 4.458 0.6990 4.9435 x 105 . 2663 0.4886 79.07 
1.4 
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Fic. 10. Curve of neutral stability for the plane Poiseuille case. 


that for the values of a for which his curve is drawn, the approximation used in de- 
riving (13.2) is no longer legitimate. Noether’s work is based upon a very good mathe- 
matical approach, which seems to promise further developments. However, in apply- 
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ing his method to particular examples, he neglected the terms in a? in the inviscid 
solutions. As is evident from previous discussions, this is bound to lead to the wrong 
conclusion that the plane Poiseuille flow is stable (as Noether actually did). The 
mathematical analysis in Langer’s work shows that the region of the c-plane for 
which c;>0 must go to zero as aR becomes infinite. This is in agreement with present 
results. Langer, however, concluded from his analysis that the motion is probably 
stable in general. This would be a natural deduction if the effect of viscosity were only 
stabilizing. The instability of the plane Poiseuille flow must therefore be attributed to 
the destabilizing effect of viscosity. This is a very significant fact and will be discussed 
in greater detail in §14. 

Pekeris’ second paper [40] is a numerical treatment of (4.1), replacing a deriva- 
tive by a ratio of two finite differences. Unfortunately, his method is not suitable for 
the purpose, because he has virtually neglected the inner friction layer. In his approxi- 
mation, he divided the half-width of the channel into (at most) four equal parts cor- 
responding to w=0, 7/16, 3/4, 15/16, 1. From the present work, we know that the 
inner friction layer occurs definitely for c<6/16. We know also from our previous in- 
vestigations that the function @ varies very rapidly in the neighborhood of the inner 
friction layer. Hence, it is not legitimate to replace ¢’ by A¢/Ay for the interval 
(0, 1/4), y being measured from the solid wall here. Also, most of the combinations 
of values (a, R) he selected do not correspond to a strong instability. These values 
are marked with crosses in Fig. 10. 

b) Stability of Blasius flow. For this case, we choose the boundary-layer thickness 
to be defined by 


y¥=6=62/V/Rz, Rt = m%/», (13.3) 
where x, § are the dimensional distances from the leading edge and the wall respec- 


tively, and %, v are the dimensional free stream velocity and the kinematic viscosity 
respectively.* With this definition, the dimensional displacement thickness is 


5, = 0.286738. (13.4) 


Such a choice has the convenience that the initial part of the velocity curve can be 
very accurately represented by 


w(y) = 2y — 3y4, (13.5) 
y being measured from the wall. Also, since the edge of the boundary layer is farther 
from the solid wall than that set by Tollmien and Schlichting, greater accuracy can 


be expected. To make it easy to compare with other results, all final values are pre- 
sented in terms of 


a, = ad}, Ry == Ré. (13.6) 


The two asymptotic branches of the a(R) curve are given by the following for- 
mulae (cf. (12.16) and (12.17)): 
1.74a, (13.7) 


4.00a;. (13.8) 


R; = 2.21(10) a; 
R; = 0.06220; , ¢ 


* Cf. Goldstein [7], vol. I, p. 135. 
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These formulae may be compared with those given by Tollmien.* The complete 
numerical result is shown in Table IV and Fig. 11. The minimum critical Reynolds 
number occurs at R,=420, agreeing fairly well with our previous estimation and the 
earlier results of Tollmien and Schlichting. 

TABLE IV. Stability of Blasius Flow. 







































































c - | a | R s | a Ry 
: 0 2.294 | 0 | x | 7214 0 oo 
0.05 2.294 0.0473 | 81.45x108 | 7214 0.0136 23 .353X105 
0.10 2.296 | 0.1040 | 4.655x10° | .7205 0.0298 1.335105 
0.15 2.311 | 0.1730 | 84555 | 7135 0.0496 24244 
0.20 2.341 | 0.2588 | 24783 ~3=| ~— 6998 0.0742 7106 
0.25 2.396 | 0.3693 | 9536 |  .6759 | 0.1059 2734 
0.30 2.481 | 0.5156 4388 | .6414 | 0.1478 1258 
0.35 2.624 | 0.7149 2358 «=| =.5897_ | (0.2050 676 
0.40 2.942 1.0778 | 1477 | = 4967, |) -0.3090 | 423 
0.411 3.21 1.2968 | 1470 | .3459 | 0.3718 | 421 
0.40 3.540 1.4264 | 1944 .3763 0.4090 | 557 
0.35 4.219 | 1.2992 | 5392 | 2893 | 0.3725 | 1546 
0.30 4.382 | 1.0055 12399 | 2733 | (0.2883 | 3555 
0.25 4.685 0.7578 34739 | = 2472, | 0.2173 9961 
0.6 
a, 
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0.4 a 
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Fic. 11. Curve of neutral stability for the Blasius case: 
—— present calculation, — — — Schlichting’s calculation. 


Earlier results. The stability of the boundary layer has been calculated by Toll- 
mien and later by Schlichting, approximating the velocity distribution by linear and 
parabolic distributions. For the evaluation of the imaginary part corresponding to the 
inviscid solutions, they used the exact profile. The calculation of Schlichting is shown 
dotted in the figure. Tollmien’s curve agrees fairly well with the present calculations, 


* Loc. cit. [73], first paper, p. 42. 
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except for a somewhat lower peak. Schlichting also calculated the amplification of the 
unstable disturbances [52], and the amplitude distribution and energy balance of the 
neutral disturbances [54]. Since the neutral curve in his calculation is inexact, it 
might be desirable to repeat some of his work if experimental results were available. 
For those calculations, the present scheme promises less numerical labor than 
Schlichting’s original work. 

14. Physical significance of the results. Prospect of further developments. Let 
us now summarize all the results which have been obtained and discuss their physical 
significance. In the first place, we may conclude that all the inertia forces controlling 
the stability of two-dimensional parallel flows can be considered in terms of the dis- 
tribution of vorticity. If the gradient of vorticity of the main flow does not vanish 
inside the fluid, then self-excited disturbances cannot exist except through the effect 
of viscosity. 

In fact, the effect of viscosity is in general destabilizing for very large Reynolds 
numbers. Thus, if a wavy disturbance of finite wave-length can exist neutrally for an 
inviscid fluid, it will be amplified through the effect of viscosity. Indeed,* if the Reynolds 
number of a flow is continually decreased, a disturbance of finite wave-length, which is 
damped at very large Reynolds numbers, becomes amplified, unless the wave-length 
is so small as to cause excessive dissipation at any Reynolds number. For still smaller 
Reynolds numbers, the damping effect becomes predominant, and we have again a 
decay of the disturbance. However, for the particular disturbance of infinite wave- 
length (essentially a steady deviation), the effect of viscosity may be said to be al- 
ways of the nature of a damping. 

The effect of viscosity is essentially one of diffusion of vorticity. It can be seen 
more clearly from the following considerations. Let us imagine a disturbance originat- 
ing from the inner friction layer where the phase velocity of the disturance is equal 
to the velocity of the main flow. During one period 27//acU of the disturbance, the 
viscous forces will propagate it side-wise through a distance of the order V/ 2rvl/acU 
=1\/27r/aRc. It is significant to compare this distance with the distance between the 
inner friction layer and the solid boundary. For if they are nearly equal, it means that 
the effect of viscosity is dominant at least from the solid surface to that layer. This 
ratio is approximately 


s = V 22/2 (14.1) 


where z is defined by (6.28). This quantity may be regarded as a measure of the ef- 
fect of viscosity. Its value is included in Tables III and IV. We notice that the value 
of s decreases from 0.7 to 0.5 as we follow the lower branch of the neutral curve of 
stability from infinite Reynolds number to the minimum critical Reynolds number. 
Then, as s decreases to zero, we are following the other branch of the neutral curve 
to infinite Reynolds numbers. Thus (see Figs. 9, 10, 11), the lower branch is essen- 
tially controlled by the effect of viscosity. The effect here is stabilizing, since an in- 
crease of Reynolds number gives instability. On the other branch, the effect of vis- 
cosity on diffusion of vorticity is predominant in comparison with the effect of dis- 
sipation. Here, an increase of Reynolds number gives stability; i.e., the effect of 
viscosity is destabilizing. This destabilizing mechanism is essentially to shift the 
phase difference between the u and v components of the disturbance. It has been ex- 


* See Fig. 9. 
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plained in some detail in Prandtl’s article [42] from the point of view of energy bal- 
ance. 

If we consider disturbances from the wall and from the inner friction layer, we 
may regard the region in between to be wholly governed by the effect of viscosity, if 
these disturbances meet after a period. Thus, it is not without significance that the 
minimum critical Reynolds number occurs for s=} approximately, which may be 
regarded as marking the passage from stabilizing effect to destabilizing effect of the 
viscous forces. 

These discussions hold both for symmetrical velocity distributions and for bound- 
ary-layer distributions. In both cases, it has been demonstrated that instability is 
essentially caused by the effect of viscosity. These velocity distributions are unstable 
whether a point of inflection occurs in the velocity profile or not. Thus, although the 
gradient of vorticity plays a part in controlling the stability of the flow, it is by no 
means the dominant factor, particularly at low Reynolds numbers. There is thus no 
reason to associate a point of inflection in the profile directly with instability. This 
removes Taylor’s objection of instability theories based on von Doenhoff’s experi- 
ments.* Even if the point of inflection in the velocity profile occurs in the leading 
part of the plate in his experiments, the flow there is definitely stable. 

There is another objection raised by Taylor against Tollmien’s work on the stabil- 
ity of the boundary layer. He questions whether the change of boundary-layer thick- 
ness should not have a drastic influence. This point can best be settled experimentally. 
So far as mathematical considerations are concerned, it seems justifiable to consider 
a boundary layer as a parallel flow; the fractional variation of thickness is very small 
over a distance of one wave-length of the disturbance, and the error incurred is only 
a few per cent. A fuller discussion of all the errors involved in the theory will be given 
in the Appendix. 

Another point should be settled by experimental investigations. Since the gen- 
eral impression had been that the plane Poiseuille flow was stable, Prandtl suggested 
that instability occurred at the entrance flow where the velocity distribution is not 
yet parabolic. The present work certainly concludes that such entrance flows are 
unstable, if they can be considered as approximately parallel. It is hard to decide 
theoretically whether a well-developed turbulence has already been reached before 
the parabolic profile is established. This presumably depends upon the conditions of 
disturbance at the inlet. The question can be best settled experimentally. 

Of the six types of problems mentioned in section 3, Part I, the three types (1), 
(2), and (5) seem to be quite settled. The present work on the boundary layer checks 
Tollmien’s result approximately, with a minimum critical Reynolds number 
R, = R6,;=420. The minimum critical Reynolds number for plane Poiseuille flow 
is found to be 10600 based upon the width of the channel. These values are at least 
not in disagreement with the existing experimental results. It would be very inter- 
esting if experiments could be carried out to check the theoretical results so far ob- 
tained. 

Since plane Couette motion is concluded to be stable while plane Poiseuille mo- 
tion is concluded to be unstable, it seems interesting to investigate a combination 
of them to find out when does the instability begin as one changes both the pressure 
gradient and the relative motion of the plates. 





* Taylor, loc. cit. [70], p. 308. 
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The stability of two-dimensional jets and wakes has never been investigated with 
the effect of viscosity included. It seems that a study of the stability of the two- 
dimensional wake might give us valuable information regarding the von Karman 
vortex street,—particularly regarding the minimum Reynolds number of its occur- 
rence and the width of the street as compared with the size of the body.* 

Transition to turbulence. The success of Taylor's theory of transition [68, 70] to 
turbulence in the boundary-layer as caused by external turbulence seems to throw 
the instability theories at a disadvantage. However, it seems that Taylor’s work 
should be regarded as only one phase of the problem, i.e., concerning cases where 
the external turbulence plays the dominant role. In fact, it is not impossible to con- 
struct an instability theory, taking account of the free turbulence outside the bound- 
ary-layer if this is the main cause of transition. The boundary condition $’+a¢=0 
at the edge of the boundary layer signifies that the disturbance there has equal mag- 
nitudes in directions parallel and perpendicular to the wall. This can be easily recon- 
ciled with the nearly isotropic turbulence in the free stream. Of course, the theory 
can only be pushed to the point where non-linear effects begin to appear. Otherwise, 
we have to deal with a very difficult mathematical problem. It is possible that the 
beginning of non-linear effect is not far from the actual point of transition. Then the 
instability theory should give useful results regarding transition, which might be ex- 
pected to check with experiment. 

APPENDIX 


In the following paragraphs, we shall describe the methods by which the numerical 
calculations are carried out. We shall then give a discussion of the numerical accuracy 
involved in the calculations. Special emphasis will be placed on the case of Blasius 
flow. 

a) Transformation of equations. The basic equations for the determination of the 
stability characteristics are given at the end of Part I. To carry out the numerical 
calculation in any particular case, we have to evaluate the functions (6.24) which 
occur in the equations (6.14), (6.15) and (6.17) through the relations (6.26). It is 
found convenient to transform (6.24) into 





ou = (1 c)(1 = a*H,)— (1 = LaMn), 
n=2 
go = Kids, — (1 — 02. a?" Non11, 
os . pa) 
giz = (1 — c)(1 — a*H2)- (ax, -> aM) + (1 — c) wi d12, 
n=2 
ore = Kidi2 + (1 — o-(1 -> aN) + (1 — c)'wd gas, 


n=l 
where the functions M,(c) and N,(c) are defined by 
M, = Hy, — H2Hy-2, n = 3, \ (2) 
Na = x. —~ KH ...1, n = a 


* This problem has been attempted by Heisenberg; see Goldstein’s book [7]. 
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The principal advantages of such transformations is to bring out the dominant terms 
in the functions di, d17 , G22 and ge. For the terms in M’s and N’s are usually very 
small (particularly for small values of c, with which we are usually concerned), while 
all the terms in the series of (6.24) are of considerable importance. This point will be 
fully discussed below. 

The calculation of (6.13) (Case 1) is still quite complicated; it is found necessary 
to take (6.15) (Case 2b) as a first approximation. Since we are not going to make any 
actual calculation for this case, we shall not go into further details with (6.13). All 
the other equations (6.14), (6.15), and (6.17) (Cases 2a, 2b, 3) can be transformed 
into the form 


_ A+trXYut+ iv) 


= —__—____ (3 
Ie) 1 + A(u + iv) 
with 7(z) defined by (6.31) and \=X(c) defined by 
wi (v1 — Yo) = — c(1 +A). (4) 


Thus, A is usually very small. The quantities u and v are real functions of @ and c, 
different for different cases. For Case 2a we have 


° / / , 
ut+tv= 1+ wiCdo22/12 


: 1 
=wic( Ki+- ) 
WC 


wic 
4-15 (1—02H2)(1—a*H2—a*Ny— --- )(Hi—a*M3—aMg— -- + +), (5) 
a 
where the second form of the right-hand side is derived by using (1). Similarly, for 
Cases 2b and 3, we have, respectively, 


ut+iv=1+ wi Ch22/¢d12 


+ w ca?(1—a?H2)(1—a*M,—a®M,— Ls )-"(N3+a?N 5+ “os ) (6) 
and 


ut iv=1+ wyic(d22+a¢22)/(d12+0¢12) 


1 
= wi e( Kit —) 
W,C 


4 wic (1—a°H2){(1—a*H2—atNy— --- )—(1—c)*%(a*®Nst+a®Nst+---)} 
a ls 4 
a? Hagih-alitd~piila <<: halls ~alie— - + ) ” 





Equation (3) contains the two real equations 
F(z) = (1+ r) {u(1 + ru) — dv?} f(1 + ru)? + (dv)?} 1, \ 
F(z) = (1 + A)of (1 + Au)? + (Av)?}-, 


for z, a, c. Thus, for each value of z, we can determine corresponding values of a 
and c. Finally, the Reynolds number is given by (cf. (6.28)) 


(8) 
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1 w; z\3 
R= sa): (9) 
we (1+A)*\ ¢ 


The actual procedure of calculation will be described presently. 

b) Procedure of calculation. The calculations required in §13 are as follows: (a) to 
find the values of @ and z corresponding to each value of c by using equations (8), 
with uw and » defined by (7) and (9); and (b) to calculate R from (9). To do this, we 
may take the following procedure. We first plot 7; against 7,; then plot the corre- 
sponding right-hand side members of (8) in a similar manner in the same diagram. 
Noting that the latter are functions of a and c only, we may plot by drawing curves 
of constant a (or constant c). The intersections of this set of curves with the (7,, 7;) 
curve give the desired results. 

This procedure is however, very laborious. A simpler method is as follows: As 
will be seen below, the imaginary parts of H’s, M’s and N’s appearing in (5) and (7) 
are very small compared with that of Ki(c), we can therefore use the approximation 


v= v(c) = — rwi{(ww"/w") for w=c. (10) 


The following steps are then taken: 
i) Calculation of aR. In this step, the auxiliary functions 
A(c), woe (c), v(c) 


are required. These are tabulated in Tables V and VI for the cases considered. Having 
calculated these functions, we can determine z and u for each value of c; aR is then 


TaBLe V. Auxiliary Functions for Calculating the Stability of the Plane Poiseuille Flow. 



































Cc we r v | wi cRl | A, | He | M; N3 
0 ia | O 0 0 0.53333 | 0.21817 | 0.06038 | 0.19340 
0.05 1.94936 | 0.01282 0.08482 0.06499 | 0.46917 | 0.20696 | 0.05047 | 0.20401 
0.10 1.89737 | 0.02633 0.18397 | 0.10187 0.41000 | 0.19516 |} 0.04183 | 0.21636 
0.15 1.84391 | 0.04061 0.30066 | 0.13015 | 0.35583 | 0.18276 | 0.03441 0.22777 
0.20 1.78885 0.05573 0.43905 0.15351 0.30667 0.16982 | 0.02813 | 0.23918 
0.25 1.73205 0.07180 | 0.60460 | 0.17356 | 0.26250 | 0.15626 | 0.02293 | 0.25000 
0.30 1.67332 | 0.08893 0.80463 0.19121 0.22333 0.14209 | 0.01872 | 0.26191 
0.35 | 1.61245 | 0.10728 1.04910 | 0.20699 | 0.18917 | 0.12729 | 0.01540 | 0.27000 
0.40} 1.54919 | 0.12701 1.35193 0.22130 | 0.16000 | 0.11182 | 0.01282 | 0.27062 
0.45 | 1.48324 | 0.14836 1.73296 | 0.23438 | 0.13583 | 0.09563 | 0.01087 | 0.26232 
0.50 1.41421 | 0.17157 2.22144 | 0.24645 | 0.11667 | 0.07875 0.00937 0.23366 
( 








determined from (9). In actual practice, a procedure of successive approximations is 
used. By taking 
( eo —_ 
WO =F), = 2, FO =F), MY = Fs) = 9 
as the initial approximation for u, z, 7,, 7:, we can obtain the successive approxima- 


tions by the formulae 


fF = LL APLC + de)? + Q0)7]™ 
Let? = FPL + rw)? + (rv)?} {CL AML + AW) J? $F AD(L + AW™) 
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TABLE VI. Auxiliary Functions for Calculating the Stability of Blasius Flow. 





c v | wf cRI H, | H; | Ms Ny 
0 0 0 | 0.60260 | 0.23513 0.07966 | 0.17615 
0.05 0.00088 0.02194 | 0.53378 | 0.22302 | 0.06975 | 0.18676 
0.10 0.00709 0.06121 | 0.46995 | 0.21212 | 0.06111 | 0.19911 
0.15 | 0.02406 0.12108 | 0.41112 | 0.19972 | 0.05369 | 0.21052 
0.20 | 0.05773 0.20469 0.35730 | 0.18678 | 0.04741 0.22193 
0.25 0.11498 0.31568 0.30848 0.17322 | 0.04221 0.23275 
0.30 0.20438 0.45835 0.26465 | 0.15905 0.03800 | 0.24466 
0.35 0.33718 0.63751 0.22583 0.14425 0.03468 0.25275 
0.40 0.52839 0.85766 0.19200 0.12878 0.03210 | 0.25337 
0.45 0.79820 1.12133 0.16318 0.11259 | 0.03015 | 0.24507 
0.50 1.17350 1 0. 0.09571 | 0.02865 0.21641 


13935 


In each approximation, Ym and 2‘ are determined graphically from 7;"’. 
ii) Calculation of a. For this purpose, the additional auxiliary functions 





wi cRi = wf eRIY KC) + - mt Hy(c), H2(c), M3(c), N3(c),--- 
WiC 

are required. These are tabulated in Tables V and VI for the cases considered. The 
methods of evaluating these functions and their accuracy will be discussed below. 
For sufficient accuracy in the final results, only the real parts of H2, M3, Ns; are re- 
quired, besides w{cR/ and H;. Having calculated these functions, we can determine 
the value of a from the real part of the equations (5) or (7). A similar method of suc 
cessive approximations may be used by writing those equations in the forms 














wic 1—a*H.—atNy— e's ’ 
a? = ——_—_ a (1 —a?H») “9 Pon=M on41/ A, 
H,(u—wj{ cRl) 1—a?P,—atPy— --- 
wic (1—a*H.—atNy— --- )—(1—c)*(a®N3g+a5Ns+ --- ) 
ga (t-~aD oa 3 Sas A 
u—wicRli (1—c)*(1—a*M,—a®M— ileal )\+a(H,;—a’?M3;—a‘M, ll 3 


An approximate value of a@ is put into the right-hand side to obtain an approximation 
of the higher order on the left-hand side. For the initial approximation, take a=0. 

c) Numerical accuracy of the calculations. The numerical accuracy of our calcula- 
tion as based upon the final equations given in section 6 are limited by several factors: 

i) by neglecting quantities of the orders e~? and (aR)~! in the reduction of the 
determinantal equations of the boundary-value problems, 

ii) by using the inviscid solutions for ¢; and @2 (error of the order (aR™')), 

iii) by the approximations of the rapidly varying solutions ¢3 and ¢,4 as discussed 
at the end of §6, 

iv) by the boundary-layer approximation used in setting up the equation of sta- 
bility (except in the cases of plane Couette and Poiseuille flows). 

Finally, certain numerical approximations have to be used in the actual evalua- 
tion of the quantities u and v in equations (11.11). We shall now discuss these factors 
one by one. 

The inaccuracy due to (i) and (ii) is negligible in all the cases considered, because 
aR is always sufficiently large. In connection with (iii), the situation is more compli- 
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cated. The first approximation of the asymptotic solution should give an error of the 
order of (aR)~'/?; while the first approximation using Hankel functions should give 
an error of the order (aR)~'/*. It might therefore be thought that the asymptotic 
method should always give a better approximation. However, this is not the case. 
For the order of accuracy of the first method is based upon a fixed value of y, while 
that of the second is based upon a fixed value of 7. Thus, if aR may be allowed to be- 
come very large while y—yo remains to be of the order of unity, the first method is 
definitely better. This is the case with the quantities de and dy. With os and 3, 
the situation is different. Here, y:—yo is always small. Except for one branch of the 
neutral curve for profiles with a flex, y1—yo goes to zero as aR becomes infinite. Be- 
cause of the smallness of y;— yo, the asymptotic solution (which fails to be accurate 
in the neighborhood of yo) never gives a good approximation. This is why the other 
method has to be used in most of the calculations, and we are limited to an accuracy 
of (aR)-'/*. We note that the curvature of the velocity distribution does not come 
into this approximation. Thus, for better accuracy, a second approximation should 
be used, the error being then reduced to the order of (aR)-*/*. However, since the 
error in the method used is only a few per cent, and an improvement in accuracy 
would not alter the general conclusions, it does not seem worth while to improve the 
accuracy in the light of general interest. Indeed, the inaccuracy due to the other 
causes (to be discussed) is also of the same order of magnitude. Another support to 
the method used is that it does agree with the asymptotic method when z is large; 
there is only a negligible difference (cf. eq. (6.31)). 

The effect of the change of the thickness of the boundary layer might be taken 
to be more serious than a mere numerical inaccuracy. Taylor regarded this as invali- 
dating the existing instability theory of the boundary layer. This question can best 
be settled experimentally. For the present, we only want to discuss its effect upon 
our boundary value problem. An approximate estimate of this effect may be obtained 
by considering the change of the thickness of the boundary layer for one wave- 
length of the disturbance. This can be easily verified* to be 4(1.72)?/a:R,. For the 
lowest value of aR; involved in the calculations of §13, this is about 6 per cent. 
Thus, the error is not large. Hence, in the physical interpretation of the results, we 
need only consider a change of Reynolds number as we pass down stream. One in- 
teresting point is the following: as the Reynolds number keeps on increasing, all 
disturbances finally become stable, if the linear theory holds throughout. Thus, the 
transition to turbulence depends upon the occurrence of the non-linear effect and 
hence must depend upon the amount of initial disturbance. 

d) Calculation of $9, ¢%, etc. We shall now discuss the method by which these 
quantities are evaluated for the calculation of u and » in the equations (5) and (7). 
A discussion of the accuracy of the present method and of Tollmien’s method of eval- 
uating these quantities will also be made. 

The original question is to evaluate the integrals H,,(c) and K,,(c) as occurring in 
(6.24). Various methods are possible for carrying out the calculation, including 
straightforward numerical integration. The method to be described is an attempt at 
a simple one. With the transformations (2), we hope to bring out the dominant terms 
of the series (6.24), and the calculations of u and v according to (7) and (9) are based 
upon the use of the transformed series (1). We make the following approximations. 


"* Cf. Goldstein [7], last column of table of p. 157. 
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i) The imaginary part v is chiefly given by that of the first term, namely, 
wi c(K,+1/wj c); this implies that the imaginary part due to H2(c), M3(c), N3(c), etc. 
is negligible. 

ii) The real part receives also little contribution from those of H2(c), M3(c), N3(c), 
etc., and hence these need be calculated only approximately. 

iii) The series are cut short; terms like Ni, Mi, --- are entirely neglected. 

Let us proceed to justify these approximations. 

The justification of (ii) and (iii) is based upon the following two facts. 

a) The quantities in the series involved decrease roughly like 1/m!, m being the 
number of integrations involved in defining a certain term. 

b) For a<1, the terms also decrease as a”. Thus, the accuracy is not very good 
for a>1, namely for low Reynolds numbers. 

But from a consultation of Tables V and VI, and the manner in which the integrals 
H2(c), M3(c), H3(c), etc., enter (5) and (7), we see that an error of ten per cent in these 
integrals will cause a negligible error in the final results. 

The justification of (i) needs more explanation. For definiteness, let us take N3(c) 
as an example. Now, 


Ve y y 
N3(c) -f dy(w — of dy(w — otf dy(w — c)~*. 
v1 


v1 ¥1 


This can be expressed as the sum of the following three integrals: 


Vo Ve 
Nac) = Kio): f dy(w — otf dy(w — c)~, 
V1 y 


Ye y ¥2 
N32(c) -f dy(w — of dy(w — otf dy(w — c)~, 
y 


Vi vo 
Ve Mv Us 

N33(c) -f dy(w — c) f dy(w — otf dy(w — c)-*. 
yo vo y 


The third integral is real, because y > yo. A further transformation of the last integra- 
tion in Nz; and N32 like 


Yo y 
f dy(w — c)-? = K,(c) -f dy(w — c)~? 
v1 


y 


Uo Vo y 
Na = {Kd}? f dy(w —c)? — Ko) f dy(w — of dy(w — c)~? 
V1 V1 v1 
Yo y 
Kio) f dy(w — of dy(w — c)? 
V1 Vo 


Yo y y 
— f dy(w — af dy(w — ot f dy(w — c)~*. 
Yo v1 


V1 


gives 


N32 


Hl 


Now, the last integral is real because y <yo. Further, it can be easily verified that 
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Vo y y v 
f dy(w — otf dy(w — c)~? -f ‘dy(w _ af dy(w — c)*. 
y v1 v1 Vo 


1 


Hence, the only term which can contribute to the imaginary part of N3(c) is 
{ Ki(c)} 2 [dy (w—c)*. Now, c is usually small so that we may put 











[are — 98 = (ny — ygr = SO 
y(w—c)?= - = — . 
Te, i hc ae 
Hence, 
Yo 1 wy ¢ 
{xo}? f dy(w —c)? = { wi cKx(c) }? — . 
v1 3 wi * 


Now we have approximately 


1 — vi. 


wi cK,(c) 


Substituting into the above expression, we obtain the imaginary part of N;(c) as 
—2w/{ cv/3w{*. This will give a contribution of approximately — (2/3) {ac(1 —c)/w{ } 2y 
to the imaginary part of v. This is negligible, because the factor preceding v is at most 
of the order of 0.02 in our calculations. Hence, it is justifiable to neglect the contribu- 
tion of N; to v. With the other terms, the approximation is even better; thus, the 
imaginary part of H2(c) is of the order of c* times that of Ki(c), and that of M3(c) 
is of the order of c* times that of Kyi(c). 

Having thus justified the approximations described above, the task is to evaluate 
K,(c), Hi(c), Ho(c), Ms(c), and N3(c) with proper degree of accuracy. For parabolic 
distribution, these integrals can be evaluated exactly; the approximation (ii) is not 
necessary. Thus, 























H,(c) = A, (11) 
Kilo) Smee {i ue} a (12) 
(c) = — —————_ + — {lo im? , 
- 2a7(1 — a?) 2a zl —a } 
He) = = =F 4 A tog (1 +0) — ogo? — — {log (1 - 09) — ix}, (13) 
3c) = — — ae — —— log a* — — jlo — a*) — tr}, 
aa ae woe Ua | loa ™ 
B2 ss 8a*A (1 a+1 
M;(c) = —- Kit { + log 
2a? 15 ta+1 a 
1 54 108 38 
+4- ~~ + =a? — at — 64ath (14) 
225 7 5 3 
1 (1 — a?*)? i+a 
RI N;3(c) = —— (1 + 3a?) — lo 
(0) = Tap S a pay 
1 t a+2x])? 
+ aif ax {(o — x*) log } on Bri ’ (15) 
16a® 0 a-— x| 
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ImN (0) T {oxi Les. oy i+a 
mN;(c) = 2B log - — lo 
, 16a° Rie 8" 3s 
2a 
ae (1 — a)(2 + 3a + 18a? + 200) (16) 
5 
where 
a? 1 8a5 
a®#=1—«, Az=at—-—+—,) B=A——-: (17) 
5 15 


These are the equations on which Table V is based, where only the real parts are 
given. For any other profile, the rule is as follows: 
i) Evaluate Ki(c) with as much accuracy as possible. Usually, it is broken up into 
two parts. Thus, 
vj 
Ki) = Kul + Kul; Kul = f “ay(w— 3, 
¥ 


Ki:(c) = fay w — c)*, (18) 


where 9; <y;<‘2. The value of y; is chosen so that Ky(c) can be calculated with suffi- 
cient accuracy by developing w as a power series of (y— yo), while K12(c) can be evalu- 
ated by developing the integrand as a power series of c/w. 

ii) Evaluate by numerical integration the quantities 


¥2 
H,(0) -{ dyw*, (19) 
V1 
v2 
H{(0) = -— 2f wdy (20) 
V1 
V2 uv 
H,(0) -{ wtdy [ wdy, (21) 


v1 v1 
V2 V2 y 

M;(0) -f widy f wtdy f wdy, (22) 
v1 = v1 
Y y y 

N;(0) -f w-tdy f widy f w~*dy. (23) 
¥1 v1 y 


iii) The integral H;(c) is then given by 
Hi(c) = Hy(0) + Hi (O)c + c?. (24) 
iv) The real part of the integrals H2(c), M3(c), Ns(c) are obtained by comparison 
with the corresponding quantities for parabolic distribution (Table V). For example, 


, 


w 2 
M;(c) — M;(0) = (=) X corresponding quantity for parabolic distribution. (25) 


The idea of the last step is essentially to approximate the given profile with a parabolic 


one. 
For the Blasius profile, (with w{ =2, y;—y, =0.4). We obtain 
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1 
Ky,(c) + —— = — 0.5615 — 0.3937¢ — 1.543c? — 1.803c* — 1.368¢* — 5.022c® 
WiC 
9 21 0.8 —¢ , 
+-:-> +2(¢ + —e+ +++) (1og——— + im}, 
8 8 c 
K 12(c) = 0.7080 + 1.3546¢ + 2.588c? + 3.860c* + 5.446c4 + 7.455c® (26) 
\ 
a 2 we 
1 
K,(c) + —— = 0.1465 + 1.2467¢ + 1.045c? + 2.039c* + 4.078c* + 2.4235 
WC 
9 21 0.8 —¢ : 
+-:: +(e + ~ ot +++) log———* + in) 
8 8 c 
In evaluating these integrals, we take 
w= 2(y— yi) — 3(y — y1)4, 0S y-—m 504, 
w=1-— {0.9—(y— y)}?, 04<y—y 5 0.9, (27) 
w= 1 09S y¥-n S11. 


For the integrals H,(0) and Hy (0), we make use of the known values of the displace- 
ment thickness and the momentum thickness 6s. 


1 
b = — (1.7208) = 0.28673, 
(28) 


1 
be +2 (1.32824) = 0.11067. 


Thus, 
H,(0) = 1 — 6; — 62 = 0.6026, H{(0) = — 2(1 — 6,) = — 1.4265. (29) 


The values of H2(0), 143(0), N3(0), as evaluated by numerical integration, are given 
in the first row of Table VI. The rest of the table is constructed by following the pro- 
cedure described above. 

We see that the method of approximation developed above is purely a numerical 
one, and the calculation can be done without excessive labor. In any case, even if the 
above method does not give satisfactory results, suitable approximations can always 
be devised for the evaluation of the necessary integrals. This is the advantage of 
using Heisenberg’s form of the inviscid solutions. In the method used by Tollmien, 
it is necessary that the profile may be approximated by linear and parabolic parts; 
otherwise, the numerical labor is excessive. It is not clear at once what is the effect 
of such an approximation on the solutions ¢(y). A more serious criticism of Tollmien’s 
method is the joining of the inviscid solutions at the point of junction of the two ap- 
proximate profiles. Mathematically speaking, such a junction presents an essential 
singularity in the coefficients of the differential equation (3.8) or (3.14). Numerically 
speaking, serious difficulty would be expected when c is equal or even only very near 
to the velocity of junction, because the inviscid solution fails at the critical layer 
where w=c. Tollmien did not publish how he overcame this difficulty.* 


* Tollmien [73], footnote, p. 37. 
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THE CYLINDRICAL ANTENNA; CURRENT AND IMPEDANCE* 


BY 
RONOLD KING anp DAVID MIDDLETON 
Cruft Laboratory and the Research Laboratory of Physics, Harvard University 


1. Introduction. The definition and the determination of the impedance of a sym- 
metrical, center-driven antenna of small, circular cross section involves three major 
problems. These are first the theoretical analysis including the formulation of bound- 
ary conditions; second the apparatus and the technique of experimental measurement; 

and third the coordination of experiment with theory. Of these 


new approach to the solution of Hallén’s integral equation in 
—+----- 0 that it replaces a function arbitrarily chosen for reasons of 
mathematical convenience in the approximate evaluation of the 
equation by a function actually fitted to the true distribution of 
current. As a consequence, new parameters are introduced to 
replace those used by Hallén (or equally those used by Gray") 
| in the successive approximations, and as would be expected the 
i - -h resulting development shows a more rapid convergence, in so 
far as this is indicated by a relatively small difference between 
first and second order solutions. 

The antenna actually analyzed is a theoretical one in the 
sense that no exact experimental analogue can be constructed. 


“2 a problems only the first is the subject of this paper; the last two 

| t are considered in detail elsewhere.! The present discussion is 

tT Ey h concerned specifically with an analytical improvement in the 

solution of the theoretical problem. 

“ oo The boundary and driving conditions in this analysis are the 

TLR, same as implied in earlier analyses,?:* and the same integral 

or 2" equation is obtained. However, the present paper introduces a 
| 
| 
| 


= -|2| 








Fig. 1. Cylindrical 
antenna with hemi- 
spherical ends. 


Its properties are summarized as follows: 
(1) The antenna is a highly conducting cylinder of small radius a extending un- 


broken from z= —/ to z= +h as shown in Fig. 1. Postulated inequalities are 
Ba K 1, a<h, (1) 


where 8 =w/c is the phase constant and c=3X10* m/sec. 
(2) The ends of the antenna at z= +h contribute nothing to the electrical prob- 


lem so that it is correct to write 
L,=0 at s=th. (2) 
(3) The antenna is center-driven by a slice generator consisting of a disk of neg- 


* Received Aug. 6, 1945. 

1 R, King and D. D. King. J. Appl. Phys. 16, 445 (1945). 

2 E. Hallén, Nova Acta, Royal Soc. Sciences, Upsala 11, 1 (1938). 
3 R. King and C. W. Harrison, Jr., Proc. I.R.E. 31, 548 (1943). 
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ligible thickness at the center, z=0, which is in all respects like any other piece of the 
antenna except that a scalar potential difference 
78 bd 
Vo = lim ($42 — 2) = O40 — -0 (3) 
z—0 
is maintained between its faces. 

(4) All other conductors and all dielectrics are sufficiently far away so that their 
individual effects are indistinguishable from the composite effect of the universe as 
a whole. If R is the distance from the center of the antenna to the nearest point on 
any other conductor or on a dielectric the following inequalities must be satisfied 


BR> 1; R> h. (4) 


The degree in which this theoretical antenna can be realized physically is sum- 
marized briefly below. Details are found elsewhere.' 45.6.7 

(a) A metal wire or rod can be constructed to satisfy completely the properties 
assumed in conjunction with (1). 

(b) If a solid cylinder with flat ends or a hollow cylinder is used (2) is not exactly 
true. A small current exists at the ends to charge the sharp edges, the end surfaces, 
or the inner surfaces of a tube near the ends. This leads to an error in h of the order 
of magnitude of a, and a consequent hidden shift in the theoretical impedance curves. 
For particular values of h near anti-resonance large errors in impedance are involved. 
A solid cylinder of length 2h along the axis with hemispherical ends as shown in Fig. 1 
is a satisfactory physical approximation that satisfies (1) and (2). 

(c) It is physically impossible to provide a slice generator. At very low frequencies 
a two-wire drive is satisfactory to approximate (3) but in this case (4) can not be 
satisfied. At high frequencies where (4) is readily satisfied a two-wire drive involves 
adjacent end surfaces and a gap in the antenna which are not taken into account in 
the theory. The effects of gap and end surfaces are compensating and may be taken 
into account roughly in comparing theoretical and experimental results by including 
a lumped capacitance in parallel with the experimentally measured impedances if the 
gap is large and a similar capacitance in parallel with the theoretical impedances if 
the gap is very narrow and the adjacent end surfaces of the antenna are very close 
together.t A vertical antenna of length h over a conducting plane, driven from a 
coaxial line, may be a good approximation of a slice generator, but the unavailability 
of an infinite, perfectly conducting plane leads to other difficulties. 

(d) The condition for the far zone (4) can not be fulfilled at low radio frequencies 
(where accurate measurements can be made easily) because it is not possible to get 
far enough away from the earth. At high frequencies where this is possible, accurate 
measurements are difficult and the dimensions of the antenna and its driving struc- 
ture become undesirably small.? 

The analysis of the theoretical antenna subject to the conditions (1) to (4) dis- 
cussed above, can be reduced to one-dimensional form involving the total current J, 
if it is assumed that the cross-sectional distribution of the density of current is inde- 


‘ L. Brillouin, Quart. Appl. Math. 1, 201 (1943). 
5 L. Brillouin, El. Communication 22, 11 (1944). 
6S. A. Schelkunoff, J. Appl. Phys. 15, 54 (1944). 
7 R. King and C. W. Harrison, Jr., J. Appl. Phys. 15, 170 (1944). 
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pendent of the axial distribution. In effect, this means that the cross-sectional dis- 
tribution and the internal impedance per unit length may be obtained from the analy- 
sis of an infinitely long cylinder. This is an excellent approximation subject to (1). 
At high frequencies the internal impedance per unit length is given by 


wull 
a (1 + 7) (5) 
2ra 


with o the conductivity in mhos per meter, yu the relative permeability, a the radius 
in meters, and II =4m X10~’ henry per meter. Subject to this assumption in addition 
to (1)—(4) the vector potential at any point z on the cylindrical surface of the antenna 
due to the axial current in the entire antenna is given by *:*-9 


ll h 

i, 2 — f Ii Roe Bidz’, (6) 
Ar h 

where 


Ri = V(z — 2’)? + @?, 


and J/ =I,(z’) is the axial current at 2’. 
The integral equation for the current, originally derived by Hallén,? is 


4a . ' 
—A, -f Tf Roe #22’ 
II =" 


10 


i z 
= == E cos Bz + 4Vo sin B\2| — vf I(s) sin B(z — sas. (8) 
Vo is the driving potential difference maintained by the slice generator at z=0; C; is 
a constant of integration which is later evaluated using (2); R.=cII =376.7 ohms 
+1207 ohms. In practice the conductivity ¢ is usually sufficiently high and therefore 
z‘ sufficiently small so that the last integral in (8) contributes negligibly to the final 
result." For simplicity it is omitted throughout the following analysis. If required 
it can be included readily at appropriate points with no change in the formulation. 
2. Expansion of the integral equation. In the absence of an exact solution of the 
integral equation (8) in closed form, an approximate solution may be obtained by 
expanding the integral on the left in a converging power series in terms of an ap- 
propriately chosen parameter. If a converging series is obtained and a sufficient num- 
ber of terms can be evaluated the choice of the parameter for expansion is unimpor- 
tant. If only a few terms in the series can be evaluated readily it is of great importance 
to select the parameter in such a way that convergence is so rapid that the sum of 
two or three terms gives a satisfactory approximation. The several parameters which 
have been used,®* including that introduced below, will be discussed critically and 
results compared in another paper. The general definition of all such parameters is 


formulated below. 


8 R. King, Electro magnetic engineering Vol. 1, McGraw-Hill Book Co., New York, 1945, p. 241. 
9S. A. Schelkunoff, Electromagnetic waves, . Van Nostrand Co., New York, 1943, pp. 140, 142 ff. 
10 Reference 3, equation (25). The complete derivation is given. 

lt R, King and F. G. Blake, Proc. I.R.E. 30, 335 (1942). 

2M. C. Gray, J. Appl. Phys. 15, 61 (1944). 
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The solution of (8) may be formulated by expressing J, in terms of a convenient 
reference current such as the input current J» and a distribution function f(z) that 


is unknown. Thus let 


I, =Iof(z); = Iz = Iof(z’), (9a) 


so that 
Tz = I,f(2')/f(z) = 1.g(2, 2’). (9b) 


The relative distribution function g(z, 2’) is defined in (9b). Now let a function V(z) 
be defined by 


A 
V(z) =f g(z, 2’) Ry 'e~ PRidz’, (10) 
= 


If the relative distribution function g(z, z’) were the actual one, it would be correct 
to write 7,V(z) for the integral on the left in (8). Whatever the form of g(z, 2’), it is 


correct to write 


4 A h 
= A; -f Ti Ree ®dz’ = I,0(z) +f [72 — I.g(z, 2’) |Rr%e~#*dz’. (11) 
~h —h 

The more nearly g(z, 2’) approximates the true distribution the smaller will be the 
difference integral on the right in (11). If g(z, 2’) can be chosen accurately enough so 
that the integral on the right in (11) is considerably smaller than the term JW (z) for 
all values of z, it is possible to treat this term as the principal part and the difference 
integral as a correction. 

If g(z, z’) were the true relative distribution function so that the difference in- 
tegral in (11) were zero, the function ¥(z) would be given by 
4n A, 


V(z) = — 
©) TI 7, 


(12) 





That is, ¥(z) would be proportional to the ratio of the vector potential on the surface 
of the antenna at a point z divided by the total axial current at z. It is clear from (6) 
that the vector potential at a point z is determined largely by the current at and 
near z, except possibly at a few points where J, is very small compared with the cur- 
rents elsewhere in the antenna. It may be assumed, therefore, that the ratio A,/J, 
is reasonably constant and predominantly real at all points along the antenna except 
at and near very small or zero values of the current. Clearly, since J,=0 at the ends 
and A, is not zero there, V(z) is infinite at z= +h. However, the product J,W¥(z) must 
remain finite and relatively small at z= +h. 

If W(z) is sensibly constant for most values of z, it must be exactly equal to V(zo) 


at some point z = 29, so chosen that Y(z) is a good approximation of ¥(z) except where 


I, is small or zero. Let 
Vv =| ¥(z0) |, (13a) 


so that 
W(z) = Vey, (13b) 


Also let a function y(z) be defined so that 
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W(z) = V+ 7(2), (14a) 


where 
y(z) = (ew — 1). (14b) 


If Y(z) is predominantly real, y(z) is a small complex correction function except at 
values of z where J, is small or zero. It is to be noted that y(z) is infinite at z= +h 
but that J,7(z) is finite and small there. 

If g(z, 2’) is not the true relative distribution function but only approximate, 
(12) is also approximate, and it is still possible to write (13a, b) and (14a, b) with 
W(z) defined as in (10). Substituting (14a) in (11) and using (11) in (8) solved for J, 
in the principal term JV, one obtains 


— jar Vas en” 6 Pon a 
I,= 2 »C; cos Bz + 4Vo sin B| 2} ; 
I \ 
— T(z) +f ag — I,2(2, 2!) [Rete Pda! ¢ (15) 
Vt h ! 


This equation is exact. Like (8) it is an integral equation in the current, but the cur- 
rent appears in the integrand of a difference integral that is small. The term J,y(2) 
is also small except near points where J, is small or vanishes, as at z= +h. 

A more useful form of (15) is obtained as follows. Let (8) be written with z=A/ in 


the form 


h 


0 = ML cd {C, cos Bh + 4Vo sin Bh} — —f Tz Ripe Phd 2! (16) 
RW (“1 2” 0 j V P* z 4XIh Ue. 
The term in z* has been omitted in (16) just as in (15). Actually (16) is exactly equiva- 
lent to (15) when this is written with z=h. In (16) 
Ry = V(h — 2’)? + a. (17) 

The desired equation is obtained by subtracting (16) from (15). It is 

— 74 , 
= - — {C,[cos Bz — cos Bh] + 4Vo[sin B| z| — sin Bh]} 


Zz 


h 


i h > bee 
~ ~ Srov(e) + f (7! — I,g(2, 2’) |Re te Ridz! - f Fi Ryxte-ttaas'h . (18) 
: ~hs —h 


This is the final exact form of the integral equation. Its principal advantage over (8) 
lies in the fact that all terms on the right involving the current are small if the rela- 
tive distribution function g(z, 2’) is correctly chosen to make the difference terms 
small. The expression (18) must be used in preference to (15) because in (18) the 
right side vanishes for all values of Bh when z= +h as required by (2), whereas the 
right side in (15) can not be made to vanish at x = +h when cos Bh=0. In this case 
the arbitrary constant C, disappears from (15). 

The integral equation (18) can be expressed as the sum of a principal current 
(I,)o consisting of the trigonometric terms and a correction current (J,), given by the 
remaining terms. The correction term (J,). can then be expanded in a power series 
in 1/¥. Thus 
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I, = (T,)o + (Tz) = (T2)o + (Ts) c + Ide + (Ts) cs + oe (19) 
Here (J,)., is obtained by substituting (Z,)o in (Zz)c; (Is), +(J2)q is obtained using 
(Zz)ot (Le) c, in (Lees etc. 
For convenience let 
F,(z) — Fa(h) = Fai Gr(z) — Ga(h) = Gaz (20a) 
where 


sin Bh (20b) 


F(z) = cos Bz; Fo(h) = cos Bh; Go(z) = sin B| |; Go(h) 


h h 
F.(s) = Fa-1.5 f g(z, 2’) Ryle #8idz’ — f Fy-1,2/Rre~ ® dz’ — Fa_s,2y(2), (21a) 
m —h 


—h 


F,,(h) 


h 
- f Fy_1,2' Rive ®*ndz’. (21b) 

—h 
The first and last terms in (21a) may be combined into F,_1,, ¥ using (10) and (14a). 
Expressions for G,(z) and G,(h) are obtained from (21a) and (21b) by writing G for 
F throughout. 

Using (19)—(21) in (18), the complete series solution for J, may be obtained. The 
constant C; may be evaluated from (16) using (19)—(21) as described in references 
2 and 3. The resulting mth order current is 


° > F,(2)/¥"- > Gah) _ > G,(z)/¥"- > F,,(h)/¥" 








12 V n= n= n= n= 
(T2)m oa tet : : : : (22) 
RY = 
Lo Fa(h)/¥" 
n=0 
This formula may be simplified using (20a, b). The result is 
(sin B(h —|2|) + D0 M,(2)/¥ 
j aV 0 n=1 
(Iz) m = , (23) 
RY 2 
cos Bh + +e F,,(h)/¥" 
n=l 
where, in particular, 
M(z) = Mi(z) + jMi (2) = F,(z) sin Bh — F,(h) sin B| z| +G,(h) cos Bz 
— G;(z) cos Bh, (24) 
M:(z) = Mi(z) + jMz (2) = F2(z) sin Bh — F2(h) sin B| z 
+ Gi(h)Fi(z) — Gi(2)Fi(A) + G2(h) cos Bz — G2(z) cos Bh, (25) 
and, as previously defined,?-*.5.6 
Fy(h) = an = Oy + jon « (26) 


With 





18 C, J. Bouwkamp, Physica 9, 609 (1942). In Bouwkamp’s paper G and F are, respectively, the F 
and G functions in this analysis. 
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B. = 8, + 78, = M,(0) = M,(0) + 7M, (0), (27) 
the impedance of the antenna is defined by 
Zo = Vo/To, (28) 
where J» is given by (23) with z=0. It is 


m 


cos Bh + >, a,/V" 


r e = ~s j c F n=1 
(Zo)m = (Ro) m oa q(x 0) m = a aaa ai e (29) 
sin Bh + >> B,/¥" 
n=] 


This is a generalization of the formula obtained by Hallén? and others*:, as shown 


later. 
3. Functions and parameters in the Hallén solution. The expressions for the cur- 
rent (23) and for the impedance (29) depend upon the constant parameter V, and 
this in turn depends upon the relative distribution function g(z, 2’). The definition of 
these quantities involves the following considerations: The relative distribution func- 
tion g(z, 2’) must be so chosen that it is a sufficiently good approximation of the 
actual distribution to make the difference integral in (11) small. Furthermore, it must 
be sufficiently simple in form that the integral (10) for ¥(z) can be evaluated and 
separated into a principal, constant, real part | W(z0)| =WV and a small correction 
term y(z) as in (14a, b). 

The choice of distribution function made by Hallén depended upon the reasonable 
albeit implicit assumption that the vector potential A, at z depends primarily upon 
the current at and near z. If contributions from all more distant elements of current 
are small, A, may be evaluated approximately by assuming the current at all points 
to be J, and neglecting retardation. This is equivalent to setting 
(30) 


§u(z, 2°) = eh), 


The subscript H will be used to designate parameters and functions in the Hallén 


analysis. With (30), (10) gives 


h dz’ h+2z ; h-—z 
WV i(z) -f = sinh™! + sinh~! ——— - (31) 
h R, a a 
Alternatively and equivalently 
Wr(z) = 20+ In (1 in ) + 6(z), (32) 
h? 
where 
2h 
G = Yes = OG) = 2In—»> (33) 
a 


as fj} [4/1 ( a ) ll / ( a ) } | 
z) =] ¥ ia are + 1 aa Tomes lees o (34) 


If the function V,(z) in (31) is plotted as a function of z/a for a range of values of 
the ratio h/a, it is found to be moderately constant for large ratios h/a except with 
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z near h. Subject to the condition a?«h?, the maximum value of Vy(z) is Q, the aver- 
age value is Q2—2+2 In 2=Q2—0.614, the value at z= +h is 32+ 1n 2. For Q>15 the 
average value or the maximum value are satisfactory approximations. In view of the 
fact that Vy(z) becomes smaller at z= +h instead of becoming infinite as it would if 
the correct distribution function were used, it is clear that the curvature of Vy(z) 
is the reverse of what it should be. Therefore, the maximum value Q is probably the 
best approximation of Vy(z) and this was Hallén’s, although not explicitly for this 
reason. Thus the Hallén analysis sets 


2h 
WV (Zo) = Wy = (= 2 In —»> (35a) 
a 


3? 

vn(2) = In (1 —_ =) + 5(z). (35b) 
1 

The Hallén expressions for the current and the impedance are given by (23) and (29) 

with Q written for VY and with appropriately modified functions F,(z), F,(h), G,(z), 

and G,(h), n>0. The functions with »=0 are independent of the choice of ¥. The 


Hallén functions are 








h e7 BR, 
Fan(z) = (Fa-1,2) 42 — (Fr-1.2)H dz’, (36a) 
= 1 
h e BRin 
Fru(h) = -f (Fn—1,2')H# dz’. (36b) 
—h Rip 


Gan(z) and G,y(h) are obtained from the above by writing G for F. These functions 
have been evaluated elsewhere*:" for n=1 and n=2. The first order distribution of 
current and the first order impedance have been calculated and represented graphi- 
cally?:!!; the second order impedance has been evaluated by Bouwkamp." The Hal- 
lén formula for the mth order current is 


yt (sin BU — |2|) + D> Man(z)/2 
n=l 








j2a4V 
( z)mH = f ; ’ (37) 
RQ = , 
cos Bh + > F,u(h)/Q" 
n=1 
is cos Bh + Zz. ann /Q" 
(Zo)mH = a oa . (38a) 
sin Bh + >> Ban /2" 
n=1 
Here 
I es : I  ; 
QnH = An + ja = Fyn h), (38b) Bau = B, + IBn “= M,n(0). (38c) 


The functions a; and §; are tabulated and represented graphically in references 2, 3, 
11; the functions a2 and 82 as calculated by Bouwkamp™ using graphical methods are 
listed in Table I and plotted in Figs. 2 and 3. 

4. Functions and parameters in the improved solution. The relative distribution 
function g(z, 2’) in (30) is the simplest and the most obvious one if an attempt is 


“ 
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TABLE I 

Bh od as B; B,' 
0 0 0 0 0 
0.2 —0.16 0.03 3.07 — 
0.4 —0.53 0.13 5.20 0.03 
0.6 —1.07 0.39 6.50 | 0.24 
0.8 —1.67 0.80 7.14 0.78 
1.0 —2.17 1.31 6.78 | 1.74 
1.2 —2.66 1.84 5.48 3.04 
1.4 —3.00 2.31 3.34 | 4.97 
1.6 3,23 2.73 0.45 7.06 
1.8 v4, 35 3.04 —3.06 | 9.33 
2.0 3.33 3.30 —7.03 | 11.81 
2.2 —3.16 3.48 —11.25 13.98 
2.4 —2.84 3.58 —16.22 | 16.08 
2.6 —2.35 3.58 —20.83 17.28 
2.8 —1.59 3.40 —24.71 17.72 
3.0 0.61 2.99 —27 .54 17.50 

2 0.50 2.37 —29.02 16.85 
3.4 1.58 1.37 —29.29 | 15.75 
3.6 2.59 0.28 — 28.30 13.84 
3.8 3.49 —0.83 —26.35 11.34 
4.0 4.33 —2.00 —23.38 8:28 
4.2 5.03 —3.09 —19.60 | 4.73 
4.4 5.41 —4.13 —15.14 0.21 
4.6 5.46 —5.04 —10.26 —4.84 
4.8 5.20 —5.67 —4.21 —10.09 
5.0 4.66 —6.08 +2.41 | —15.10 


made to solve the original integral equation (8) as was done by Hallén. On the other 
hand, if the formal solution is carried through to obtain (23) without previously 
selecting g(z, 2’) as has been done in the present analysis, it is perfectly clear that the 
leading term in the distribution of current for any value of the distribution function 
must be of the form 
I, = Kfi(z) (39a) 
with 
fi(z) = sin B(h — | 2| ). (39b) 
K is an amplitude factor independent of z. Accordingly, an approximate relative dis- 
tribution function is 


sin B(h — | 2’|) o fale’) ' (40) 


sin B(h — | z| ) fiz) 

This function is known to be a very much better approximation of the actual current 
then the function assumed by Hallén, gu(z, 2’) =e#*. The function fi(z) =sin B(h — | | ) 
actually is proportional to the principal part of the current; the function e*: is not. 
Using (40) and (10) we obtain 
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Fic. 2. The parameters a} and a,’ as a function of Bh. 





Fic. 3. The parameters 8, and a as a function of fh. 
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h 
Vxi(z) -{ gxi(2, 2) Rye PRid2’, (41) 
h 


This function involves the factor f(z) =sin B(h—|s|) in the denominator of the 
integrand. Since this is not a function of z’ it is a constant in the integration. There- 
fore, it is convenient to introduce the function 





¥i(z) = f * fle)Rrte- Rida! (42a) 
so that ° 
fede o (42b) 
f(z) 
The function ¥;(z) can be written in the form 
V,(z) = C(z) sin Bh — S(z) cos Bh, (43) 
where 
C(z) ={ vn Bz’ Rye #idz’, (44) 
- 
S(z) ={ sin B| 2’ | Ry'e— #¥1dz’, (45) 
—h 


These integrals are evaluated in the Appendix both in general and in a simpler ap- 
proximate form. The latter is a good approximation if, as assumed throughout this 
analysis, h*>>a?. Curves for C(z) and S(z) as calculated using the simpler forms which 
apply in this analysis are given in Figs. 4-7, 20-23 for Bh=a/2 and mw and for 2 
=2 In (2h/a) =10 and 20. It is to be noted that 


¥i(z) =C(z); Bh=-x/2, (46) 
¥i(z) = S(z); Bh =r. (47) 


It follows that the plots of C(z) with 84=7/2 are also plots of ¥(z); these are given 
in Figs. 4 and 5 for 2=10 and 20. Similarly plots of S(z) with Bh=7 are also plots 
of ¥,(z); these are given in Figs. 6 and 7 for Q=10 and 20. The function y;(z) is seen 
to have a very small imaginary part so that it and Vxi(z) =¥(z)/fi(z) are predomi- 
nantly real, in confirmation of the assumption made in conjunction with (14). Ac- 
cordingly, the parameter V = | W(zo)| defined in (13a) may be chosen to be 





Vxi =|¥xi(0)| =|¥(0)|; BA = 2/2; (48) 
Wai =| Uei(h — d/4)| =|vi(h —2/4)|; Bh =o. (49) 
The function (a) | 
, | alz) | ‘ 
Vv c1(z = - 0 
| Vxi(z) | file) (50) 


is plotted in Figs. 4-7. For Bh=2/2 and both for 2=10 and 20 it is seen to be quite 
constant over the entire length of the antenna except near the ends where it becomes 
infinite, as it should. For Bh=7 the function becomes infinite not only at the ends 
but also at the center. The infinity at the center is a result of approximating the 
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current by the distribution function f(z) =sin B(h— | z| ). With Bh=a and z=0, this 
vanishes so that Vx;(z) necessarily diverges. Unlike the infinity at the ends, the in- 
finity at z=0 is due to the fact that f;(z) and hence gx(z, 2’) are approximate and not 
exact distribution functions. Actually, the current does not vanish at z=0; it merely is 
small so that Y(z) does not become infinite. The fact that J, is small at and near 
z=0 does not mean that V(z) necessarily becomes very large. V(z) is by definition 
proportional to the ratio A,/J,, and A, is determined largely by the current at z. 
Hence A,» is small if J,» is, and the ratio may remain moderately constant. Fur- 
theremore, since A, at z=h—}/4 is determined principally by the-large (near maxi- 
mum) currents at and near z=h—X/4, it is affected only very slightly by a small 
current at z=0. Therefore A, at z=h—X/4 and Vxi(4—X/4) will not be sensibly dif- 
ferent if a fictitious zero current is assumed at z=0 or an actual small current. Ac- 
cordingly the function | Vxi(h—d/4| is a good approximation of Vx1(z) for the actual 
current everywhere (including z=0) except near the ends, z= +h. 

Although the qualitative argument to show that Vx:(z) is sensibly constant and 
approximately equal to | Wi(h—d/4)| for all values of z except the ends is sound, 
it can be verified directly using Hallén’s first order distribution. It has been shown’ 
that a very satisfactory approximation of the Hallén first order current is given by 





cos Bz — cos Bh 
I, = Ig’ ( 

1 — cos Bh 
where Jj’ is the component of current at z=0 in phase with the driving potential 


difference and J,’ is the maximum value of the component of current in phase quad- 
rature with the driving potential difference. Im occurs at z=h—)/4. With 


) +t sin B(h — | z| ); = S Bh < 2r, (51) 


k= I6'/Tel; |k| <1, (52) 


it is possible to write (51) in the form 


cos Bz — cos Bh ; 
I, = jl {sin B(h —|z|) — ir( p = jl fo(z). (53) 
1 — cos Bh 





With this approximate current, an appropriate distribution function g(z, 2’) is defined 


by 


|B sin B(h — | 2’ | ) — jk(cos Bz’ — cos Bh)/(1 — cos Bh) 











g2(2, 2) = I. sin B(h— | z|) — jk(cos Bz — cos Bh)/(1 — cos Bh) 
J2(2’) ul 
= : — S Bh < 2r. - 
F(a) 7 = Bh < 2 (54) 


The ratio factor k is negative and small compared with unity. It is plotted in Fig. 8 
as a function of Bh from the data of Figs. 9-11 in reference 3. Only values of Bh near 
7 are used because for Bh not near integral multiples of 7 the distribution (cos Bz 
—cos Bh) does not differ greatly from sin B(h— | s| ). At Bh = 7/2 they are identical. 
Using the notation (42a, b), 
h 


¥2.(z) = fo(2’) Rr te #*1d2’, (55a) 
al 


so that 
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Vx2(z) = ‘ (S5b) 


The function W2(z) can be written 





ee _ [C(z) — E(z) cos Bh 
Yo(z) = C(z) sin Bh — S(z) cos Bh — i] | (56) 
1 — cos Bh 
where C(z) and S(z) are defined in (44) and (45), and E(z) is given by 
h 
E(z) ={ Rye #Ridz’ (57) 
= 





Fic. 8. The quantity k=I}’ /I,/ as a function of 8h near anti-resonance. 


This function is evaluated in the Appendix both in general and in a simpler ap- 
proximate form valid when it is possible to write a?<h?, as in the present analysis. 
E(z) is plotted in Figs. 24 and 25 for Bh=7/2 and mw and with 2=10 and 20. The 
function 2(z) is necessarily predominantly real because it is known that the first 
two terms in (56)—these are identically y,(z)—are predominantly real and that k is 
small. The functions | Po(z)| and | Vx2(z)| = | P2(z)| / fe(z)| are shown in Figs. 9 and 
10 for Bh=7 and Q=10 and 20. It is seen that | Vxe(z)| does not become infinite at 
z=0, and is reasonably constant and equal to | Wxe(h —d/4)| = | ¥2(h—2/4)| for all 
values of z except at the ends where it becomes infinite, as it should. Comparison of 
Figs. 9 and 10 with 6 and 7 shown that | Po(h —d/4)| differs only slightly from 
¥i(h—\/4)|. The difference is greater for the smaller value of 2. It follows that 
¥i(h—/4)| is a satisfactory parameter even for Bh =7. If desired | W2(h—d/4) | may 
be used especially for small values of 2, but the difference is not over about 3% for 
Q= 10. 
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The conclusions of the above analysis may be generalized and summarized as 

follows: 1. The relative distribution function 

sin B(h — | 2’| ) 

sin B(h — | 2 | ) 


is a good approximation for all values of h. 2. Suitable parameters for expansion are 





g1x(z, 2’) a 





Fic. 11. The parameter y, Eqs. (48) and (49), as a function of Bh, 
for 2=10, 15, 20. 


| Wx1(0)| =|yu(0)| for BhSx/2; | Wxi(h—2/4)| =|yr(h—2/4)| for Bh2 1/2. The fol- 

lowing notation will be used from here on 

Y { ¥xi(0) | = | ¥1(0) |; Bh < x/2 
| Wxi(h — /4)| =|vlh — 2/4); BAS 2/2. 

Since Vx:(z) has such a small imaginary part and is so well represented by (58) ex- 

cept at the ends of the antenna, the correction function y(z) in (14a) is sufficiently 


small to be neglected. 
The parameter y as defined in (58) is plotted as a function of Bh for 2=10, 15, 20 


v (58) 
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in Fig. 11. For Q=10 the curve in solid line for Bh>~2/2 is | Y2(h—d/4)|, the curve 
in broken line is |¥i(h—X/4)|. The curves for 2=15, and 20 with Bh>w/2 are 
| ¥i(h—-2/4)|. 

5. Distribution of current. The distribution of current is given by (23), the im- 
pedance by (29) with W(=yW) obtained from (58) or Fig. 11 for the value of h/a in 
question. The functions F(z), Fo(h), Go(z), and Go(#) are unchanged; they are de- 
fined by (20b). Upon substituting (40) in (10) and using (10) in the form (14a) in 
(21a) with »=1, this becomes 

h 


Fix(z) = FosWxi(z) — Fo. Ry'e~ #®1dz’, (59) 
—h 


Upon comparing (59) with (36a) written for »=1, it follows that 


Fix(z) = Fy n(2) +} (y = 2) (Foz) #. (60) 
Since V(z) is not involved in (21b), 
Fix(h) = Fin(h). (61) 
Upon substituting (60) in (21a) with y(z) =0 and »=2, this becomes 
A 
Fer(z) = (Fi2) KV xi(2) - (Fie) KRy'e~ P*idz’. (62) 
ol 


Using (60) and (61) in (62), this gives 
h h 
Fox(z) = (Fis) xw -f (Fi) nH Rye #Fidz’ — (yy — » f Fo, Rye *#* dz’. (63) 
wth —h 
Upon comparing (63) with (36a), this time written with n=2, and using (60) and 


(61) as well as (36a) written with n=1, we see that 
Fox(z) = Pis)ad + WW — 2) (Foe) + Fan(2) — (Fis) a2 
-~ = 2) [(Fos) av — Fin(z) — W- 2) (Fos) x]. 
When terms are collected there results, 


Fox(z) = Fen(z) + W — 2)(Pis)a + (WY — 2)Fin(s) + (W — 2)?Fos) x. (64) 


Using (21b) with (60) and (61), we find - 
Fex(h) = Fen(h) + (W — 2)Fin(h). (65) 
Subtracting (65) from (64), we have 
(Fos)x = (Fas)a + 2 — 2)(Fis)a + (W — 2)?(Pos) a. (66) 


Repetition of the above procedure to evaluate (F,.)x leads to: 


n 


! 
Pde 0 Sh OPnadet: -0Ee (67) 
imo (" — jit! 


Expressions for Gex(z), Gex(h), (Ges)x, and (G,z)x are obtained from (64)—(67) by 
writing G for F. 
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If the functions F and G are combined to form the functions M as defined in (24) 
and (25) for the first two terms, the result is 


Mix(z) = Min(s) + (YW — Q) sin B(k — | 2] ), (68) 
Mox(z) = Mon(z) + 2(¥ — 2)Min(z) + (¥ — b)? sin B(h — | 2] ). (69) 
In general, 
2 n! 
M,x(z) = >» — (y — b)*M,_:,n(2); n= 0, (70) 





i=0 (n — 4) !i! 
where it is understood that 
Mou(z) = sin B(h —|2|). (71) 
Similarly, 
n—! (n — 1)! 
F,x(h) = =. (y — Q) Fi n(h); n 


ino (2 — i — 1)!2! 


IV 


(72) 





Upon substituting (70)—(72) in (23) the general expression for the mth order 
current becomes 


m 


¢{ (Di)m sin B(A — | |) + DS (Dass) mM nu(z)/¥" 








iiewtoe = (73) 
Rwy = 
cos Bh + >> (Dn) m—1F nn(h)/y* 
where with ii 
x=1- had (74) 
v 
the D’s have the following significance: 
Order m = 0 | 1 | 2 | 3 | 4 | 
(D)m={ 1 | +2 +22 | ta* | txt | +-- 
(Ds)m = | 1 | + 2x + 3x? + 4x3 +-- 
(Ds)m = | | } 1 +3e | +6x* | +--- (75) 
(Dim = | | tf | +42 | +-- 
(Ds)m = | | 1 ios 





It is interesting and significant to note than when the series in (75) are summed 
for an infinite number of terms, i.e., m—>o, then 


qr" 1 1 v\" 
D, =- =( ) =: = (~) :  < 4; n= 1. (76) 
dx*™'\1— x (1 — x)" Q 


With these values of (D,)» and an infinite number of terms, (73) is identical with the 
expression (37) obtained by Hallén. Furthermore, if a—0, Y— for all values of Bh, 
so that (73) approaches (37) as the radius a approaches zero. 

It is important to note that if a finite number of terms is used in (73), all terms 
belonging to a given order m of solution must be retained and no others. That is, if 
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Fic. 14. First order current for Bh=2/2 in units of Vj/602Dy. 
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Fic. 15. First order current for Bh =3x/4 in units of V>/602Duz. 
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an mth order solution is evaluated, only terms contributed by M,x(z) and F,x(h) 


with n=0, 1, 2,--+-+, m are used. It is readily verified that this is equivalent to 
writing 
ovat s(1-2) (0-8) 
| Fes pce 6 (ck ae 
¥ v 
2 
(Dis = 14 2(1 “ -), (77a) 
¥ 
(Dz) = :.. 











04 06 os 





02 
Fig. 16. First order current for Bh=-= in units of V;/600Dyz. 


The expressions (D;)2 and (Dz)2 are plotted as functions of Bh for 2=10, 15, 20, in 
Fig. 12. Similarly it is correct to set 


(D,)o = ., 

2 
(Di)1 = 1+ (1 - =). (77b) 
(D2), = 1. 


The function (D;); is shown in Fig. 13 for Q=10, 15, 20. 

The first order distribution of current as calculated on the one hand from the 
Hallén formula (38) in reference 3, and on the other hand, from (73), are shown in 
Figs. 14-16 for Bh=7/2, 32/4, 7, and with Q=10, 20. The function f” and f’ when 


(73) is written in the form 


By a 2xVo 4 ae (78) 
(1 = Sop. Y if’) 
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are plotted in the figures. Numerical values of Dy are given in reference 3 where D 
is written instead of Dy. 

Apart from the change in the input current which is discussed below in terms of 
the impedance, the general shape of the two sets of curves is much the same. The 
new, more exact theory leads to a distribution with somewhat greater relative ampli- 
tudes nearer the outer parts of the antenna, and with a somewhat larger component 
in phase with the driving potential difference. For 8h = 7 the first order solution of the 
new theory is the same as the first order solution of Hallén’s theory if y is written for 
Q. Since W is less than Q, this means the new first order distribution is the same as 
Hallén’s first order distribution for an antenna of greater radius, but only for Bh =r. 

6. The impedance. The formula for the impedance according to the new, more 
exact theory is 

m 
cos Bh + >> (Dn) m—10n/¥" 
— Je n=1 


_* (79) 
2r : m 
(Dy) m sin Bh + >> (Das1) mBn/¥" 


n=l 





where a, and £, are defined in (38b, c); a; is tabulated in reference 11; a2 in Table I. 
Curves for (Ro)m and (Xo)m as calculated from (79) are given in Figs. 17-19. Both 
second and first order solutions are shown for 2=10, 15, 20. Thesé are calculated 


from 

A, + 7A32!| 

it, ii 
By, + 7B: 


7 (tan-1A42/A1—tan-1B2/B}) | (80) 





where for the second order solution 


Ay = (Dy)1a1 /¥ + (D2)102 /V’, 
Ag == [cos Bh + (D;)01/ a (Dz) 103/¥'], 
By = (Dy)2 sin Bh + (D2)2Bi/v + Ba/¥ , 
II 1, 2 
By = (D2)2Bi /¥ + B2/¥ , (81) 
with (D,)2 and (Dz): given by (77a) and (D;):, (Dz): by (77b). 


For the first order solution 
A, = (D1)oat » 
A, = — [¥ cos Bh — (Ds)oas], 
B; = ¥(D;): sin Bh + Bi, 
B, = Bi, (82) 


HT 


where (D;)o and (D;); are given by (77b). 

The impedance calculated according to the new, more exact theory differ con- 
siderably in some details but not in major outline from that obtained from the 
Hallén theory as calculated by King and Blake,’! King and Harrison,’ and Bouw- 
kamp.™ In general, resistances at antiresonance are smaller and occur at smaller 
values of Bh; resistances at resonances are greater and likewise occur at smaller 
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values of Bh. These differences are most significant for large values of the radius of 
the antenna. A critical discussion of impedance calculated from the theory here 
presented and comparison with experiment and with the theories of Hallén,?.* Gray,” 
Schelkunoff,*:' and others is reserved for a sequel to this paper. Accuracy of the re- 
sults and convergence of the series involved also will be discussed therein. 





Fic. 19. The input resistance and reactance of a very thin, cylindrical, 
center-driven antenna, 2=20, h/a=1.1 X10‘, for the first and second order 


theories. 


Using a specially constructed coaxial line and driving conditions that approxi- 
mate as closely as possible the idealized slice generator D. D. King has measured 
the impedance of cylindrical antennas with hemispherical ends. A complete descrip- 
tion of the apparatus, of the method, and of the results will be contained in a doctoral 
dissertation and in a paper to be published in another journal. A cross-section of the 
results involving all of the critical values for 2=10 are shown below together with 
the corresponding theoretical results of the theory outlined above. The agreement is 
seen to be good for ail quantities in the second order theory, only approximate in the 
first order theory. 
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TABLE II 
Anti + | 
; 4 X min | 
resonant | | - [m — |¢—Bhanti-res. Resonant = Bhres. |Ro at Bh= id X,at ph= - 
re Xie | Ro 2 2 2 








Experimental | 
Results by 800 | 1.95 | .60 2.5 | SY 85 47 
D. D. King | 








Theoretical | 
Results 860 1.80 61 71.0 .094 88 42.5 
Second Order | 








Theoretical | | 
Results 840 1.27 .39 64.8 -065 67 
First Order 
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APPENDIX: INTEGRAL FUNCTIONS 
The Functions C(z) and S(z). 


h h 
C(z) = J R7'e—#*®1 cos Bz'dz’ = f (Ryte~ #81 + Ry 'e—22) cos Bz'dz’, (1) 
h 0 
h h 
S(z) = f Re'e—#* sin B| 2’ | dz’ = f (Ry te #*® + Ryte—#R2) sin Bz’dz’, (2) 
will 0 
where 
Ri = V(z— 2)? + @ (3a) R. = V(z + 2)? + a’. (3b) 
These integrals can be written in the form 
Cz) =3{hth+1+ i), (4) 
ae j 
S(z) «—h-h+h-he lh hth hl, (S) 
where 
h h 
I; -f Rye BRi-2") da! - cm f Rye BR?) dz’ (6) 
0 0 
h h 
I, -f Ry te 8 Rirt2" dz! = e— ibe J Rye Bist 20g! (7) 
0 0 
A A 
I; -f Re 'e~ 8 (R22) dz! = ¢e~ ibs Ry 'e7 #8 (Bs 2-2") dz’ (8) 
0 0 
h A 
I, = f Ry le— iB (Rete) dz! = eibe f Ry te~ 8 Ratete’) de’, (9) 
0 0 


The four integrals (6)—(9) can all be reduced to the form 


% 
f v—e—*dv (10) 
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by making the changes in the variable and in the upper and lower limits listed in 


Table III. 














TABLE III 
dv 
Integral ae vo for 2’=0 | » for 2’=h 
| Ss z v 
| 8(Ri+2—-2')=j8[V G2)" +0? (@'—2)] le] =- — | i8(Rots) | 58(Ru—m) 
Ri 1 
I, | jB(Ri-— +2’) =78[V/< (2 (2’—2z)?+-a?+(2’— 2) aes +1] == jB(Ro—2) j8(Riut+m) 
I; j8(R2—2—2') =j8[/(@’+2)*+a*—(2’+2) ] of Ett | =—— | jB(Ro—z) | j8(Ren—w) 
pee » 
I, | §8(Ret2+2') =j81V/ +2) +47+ (2'+2) | ial, +1] =~ | i(Rot2) | 58(Ratus) 
2 2 | 
Ro=v z?-+a?: u2=h+z; u=h—z; Ra=-v u2+a?; Ru=V/ui+a? 


In terms of exponential, sine, and cosine integrals, 


u 


f © du = Ei (— jb) — Ei(— je) = Ci (6) — Ci (a) —f Si(b) +7 Sia). (11) 


With (11) the several integrals (6)—(9) may be expressed as follows in terms of the 

exponential integral and the sine and cosine integrals 

T, = — e%:{Ei [— j8(Ru — m)] — Ei [— j8(Ro +2)]} (12a) 
= — ¢6:{ Ci B(Ry, — wu) — CiB(Ro + 2) — 7 SiB(Rin — u) + 7 SiB(Ro + z)}. (12b) 


I. = e~#*{Ei [— j8(Ris + m)] — Ei [— j8(Ro — 2)]} (13a) 
= ¢~#8?{ CiB(Ri, + m1) — CiB(Ro — 2) — fj SiB(Rin + m1) +7 SiB(Ro — 2)$. (13b) 
Iz = e~#*{ Ei [— j8(Rex — u2)] — Ei [— j8(Ro — 2)]} (14a) 
= ¢-i8:{CiB(Ro, — u2) — CiB(Ro — 2) — 7 SiB(Ren — u2) + 7 SiB(Ro — 2)}. (14b) 
I, = e%*{ Ei [— 78(Ron + u2)] — Ei [— j8(Ro + 2)]} (15a) 


eit Ci B(Ro, + us) — CiB(Ro + 2) — j SiB(Ren + m2) +7 SiB(Ro+2)}. (15b) 


Il 


Upon combining the several integrals according to (4) and (5), 
C(z) = fei*{ Ei [— j8(Rea + u2)] — Ei [— j8(Rius — m)]} 
+ 4¢-#2{ Ei [— 78(Rin + m1)] — Ei [— j8(Rea — u2)]}. (16a) 
C(z) = 4e%*{ CiB(Ren + u2) — CiB(Rin — 1) — 7 SiB( Rea + m2) + 7 SiB(Rua — uy) } 


+ 1¢-i8:{ Ci B(Ri, + 1) — CiB(Rea — U2) — 7 SiB(Rin + m1) 
+ 7 SiB(Ren— u2)}. — (16b) 
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S(z) = = om Ei [— j8(Rea + u2)] + Ei [— j8(Ris — u1)] — 2 Ei [— j8(Ro + 2)]} 


+ 2 e~ #8! Fi [78(Rir + u)] + Ei [—j8(Rer — u2)]— 2 Ei [—j8(Ro — z)]}.(17a) 


< 


S(z) = : e82(Ci B(Ro, + U2) + Ci B(Riu, “—~ 1) = j Si B(Re» + U2) 
— 7 SiB(Rin — uy) — 2 Ci B(Ro + 2) + 72 Si B(Ro + 2)] 


+ = e-#*[Ci B(Rin + 1) + Ci B(Ron — 2) — 7 Si B(Rin + mi) 


— 7 SiB(Ren — us) — 2 Ci B(Ro — 2) + j2 Si B(Ro — 2)). (17b) 
In trigonometric form 
C(z) = 4 cos Bz[Ci B(Rea + ue) + Ci B(Rin + 1) — Ci B(Ren — ue) 
— CiB(Ri, — m4) — 7 SiB( Ren + u2) — 7 SiB(Rin + u1) +7 SiB(Rea — U2) 
+j Si B(Rirn = u;| 
+ 3 sin Bz[Si B(Rex + 2) — Si B(Rin + m1) + Si B(Rex — 2) 
— Si B(Rix — uw) + 7 Ci B(Rea + ue) — 7 Ci B(Rin + m1) + 7 Ci B(Rern — ue) 
—jCiB(Ru—m)]. (18) 
S(z) = 4 cos Bz[Si B(Rex, + we) + Si B(Ria + 1) + Si B(Rea — m2) 
+ Si B(Rin — wu) — 2 Si B(Ro + 2) — 2 SiB(Ro — 2) +7 Ci B(Rea + m2) 
+ 7 Ci B(Rin + m1) + 7 Ci B(Ren + ue) +7 Ci B(Rin — 1) 
— j2 Ci B(Ro + 2) — j2 CiB(Ro — 3)] 
— } sin Bz[Ci B(Re, + u2) — Ci B(Rin + m1) — Ci B(Rer — we) 
+ Ci B(Rin — uy) — 2 Ci B(Ro + 2) + 2 Ci B(Ro — 2) 
— 7 SiB(Roa + u2) + 7 SiB(Rin + u1) + 7 Si BCR — U2) 
— j Si B(Rin — m1) + j2 Si B(Ro + 2) — 72 Si B(Ro — 2)]. (19) 
With 
R=VP+a; wm=ht+s; m=h—-s; Rn= VFO; Rn = Vat or 


These are exact expressions for the integrals (1) and (2). They are valid for all values 
of the argument z and the parameters h and a. 


A ¢@— ibR1 
E(s) = f dz’. (20) 
-» 


The integral E(z). 





Let the variable be changed by setting 


BR: = BV (z — 2’)? +a? = VU? + V?, (21) 





where 
U = B(z — 2’); V = Ba. (22) 
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The integral then becomes 


BR» 
E(z) -f (U2 + V2)-%¢-4 07+?) 1/2g 
F 


Rip 


BRoh 
~ f (U2 + V2)-"/2 cos (U2 + V2)""dU 
8 


Rip 
BR», 
_ if (U? + V?)-!/? sin ]U? + V?)'/*7dU (23) 
BRip 
with 
Ran = V(h +2)? + a; Ruy = V(h — 2)? + a. (24) 


Let the following symbols be introduced: 


Cuv x -f (U? + V?)-!/2 cos (U2 + V?)!/*dU, (25) 
0 
Suv x -f (U2? + V?)-!/2 sin (U? + V?)!/4dU. (26) 
0 
These functions satisfy the conditions Cuv (—x)=—Cuv x, Suv (—x)=—Suv x. 
In terms of the notation (25) and (26), the integral (20) becomes 
E(z) = Cuv BR, — Cuv BRi, — 7 Suv BR2, + 7 Suv BRin. (27) 


This is an exact expression for the integral (20). 


Approximate Expressions for C(z), S(z), and E(z). 


compared with h, useful approximate expressions for the integrals C(z), S(z), and 
E(z) may be derived as follows. Expanding the integral cosine using 


Cix=C+Inx—- cia= f wi — cos u)du (28) 
0 

where C is Euler’s constant, one obtains 
| Ci B(Ren — us) = C + In B(Ren — 2) — Ci B(Rea — m2), (29) 
Ci (Rin — m1) = C + In B(Rix — 1) — Ci B(Riu — 1), (30) 
Ci B(Ro — z) = C + In B(Ro — 2) — Ci B(Ro — 2). (31) 

However, 

Roy — U2 = V(h + 2)? + a? — (h +2), (32) 
Ruy — 1 = V(h — 2)? + a? — (h — 2), (33) 
Ro —-2s=V2+a—z, (34) 


are all of order of magnitude a, so that the arguments B( Re, — U2) and 6(Ri,— 4) are 
of magnitude fa. It follows that since the arguments of Ci B(Re,—u2), Ci B(Rins— 41) 
and Ci 8(Ro—z) are small these functions may be expanded in series. The leading 
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term in each case is of the order of magnitude 38%a? so that the Ci terms in (29)—(31) 
are negligible compared with the logarithm. Hence, 


Ci B(Ro, — ue) = C + In B(Ron — 2), (35) 
Ci B(Ri, — m1) = C+ In B(Rua — 1), (36) 
Ci B(Ro — z) = C + In B(Ro — 2). (37) 


Since the functions Si B(Ry~—u) and Si B(R2,—u2) are of order of magnitude Ba, 
they are negligible compared with Si B(Re+-u2) and Si B(Ri+u1) except very near 
the ends z= +h. If Ci B(Ra+u2), Ci B(Ru+u,), and Ci 8(Ro+2) are expanded using 


(28) and the relation 























u + (u? + a?)1/? u 
In = sinh~! —, (38) 
a a 
and the approximations, 
Ci B(Re» + U2) = Ci 2Bue2 => Ci 2B(h te Z), (39) 
CiB(Ris + ux) = Ci 26u, = Ci 28(h — 2), (40) 
CiB(Ro + 2) = Ci 26:, (41) 
(18) and (19) reduce to the following approximate forms: 
C(z) = — 4cosBz[Ci2B(h + z) + Ci2B(k — 2) + 7 Si2B(h + 2) +7 Si 2B(h — z)] 
+ 4sin Bz[Si2B(h + 2) — Si2B(h — z) — j Ci2B(h +2) +7 Ci2B(h — 2)] 
h+2 h—z 
+ cos Bz sink 1. +.sinh-!— |. (42) 
a a 
S(z) + 3 cos Bz[Si 28(h + z) + Si 28(k — 2) — 2 Si 28| z| 
— 7 Ci 26(h + 2) — j Ci 2B(h — 2) + 27 Ci 26] 
+ 2 sin Bz[Ci 28(h + 2) — Ci 28(hk — z) + 7 Si 28(A + 2) 
— 7 Si 28(h — 2) — 27 Si 282] — sin B| z| Ci 262 
: ‘ey ee h+z : h-z 
+ sin B| z|| sinh! + sinh 
a a 
, h+|z| z | 
— 2sin 8} z} sin A sinh! | -|. (43) 
a a 


The last factor in (43) is written in the expanded form shown in order to contain 
Wz(z) =sinh—(h+z/a)+sinh—(h—z/a). The remaining two terms may be written in 
the following approximate form if desired, 


h+|s| 2 | 
- sink —— = sinh™ oe 


a a 


' h+|2|+ [(A+2)?+ a*]'” 
= — ln —— — —— 


| | + (cg? + a?) 


= In( || ), (44) 
h+|z| 


The function (42) is shown graphically in Figs. 20 and 21 for Bh=7; in Figs. 4 and 
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5 for Bh=7/2. The function (43) is in Figs. 22 and 23 for Bh=7/2; in Figs. 6 and 
7 for Bh=r7. 

An approximate expression for E(z) is obtained from (23) by adding and sub- 
tracting 


BReh 
f (U2 + V2y-1aU, 
B 


Rin 


BReh 
E(z) = -f (U2 + V2)-2[1 — cos (U2 + V2)!/2]dU 
B 


Rih 


BReh 
po if (U? + V2)-1/2 sin (U? + V2)!/2dU 
B 


Rip 
BReh 
v (U2? + V*)-"/*7dU. (45) 
BRiy, 
If the small quantity V =8a is neglected in the first two integrals in (45) these re- 
main everywhere finite and vanish at U=0. This is not true in the last integral in 
which V plays an important part. If V is neglected in the first two integrals both in 
the integrand and in the limits, but retained in the last integral, the following ap- 
proximate expression is obtained: 


B(zt+h) B(2zt+h) 
E(z) = — | U|-(1 — cos U)dU — if U- sin UdU 
“ B(2z—h) B(z—h) 
BR 
+ (U? + V?)/*dU. (46) 
BRiy 


Because the magnitude of U and not U itself appears in the denominator of the first 
integral, this must be evaluated in two steps for the ranges 2’ >z and 2’ <z. It is not 
necessary to write | U| in the second integral because the integrand does not change 
sign as 2’ passes through z. Hence with U =8(z’—z) = —U the first integral in (46) 
becomes 


B(z2+h) 0 ae pe ere 
-f U-"(1 — cos U)dU -f U-(1 — cos U)dU 
0 


B(2z—h) 
B(h+z) 
= -f U-(1 — cos U)dU 
0 
B(h—<) oes par 
-f U-"(1 — cos U)dU. (47) 
0 


Using (28) in (47), the sine integral in the second integral in (46), and evaluating the 
third integral directly, we write (46) in the form: 


E(z) = — CiB(h + 2) — CiB(h — 2) — 7 SiB(A +2) — 7 SiB(h — 2) 
h+z2 h— 
+ sinh (- ~*) + sinh ( *). (48) 


a a 





Use has been made of the fact that Si B(z —h) = —Si B(h—z). The function E(z) in 
(48) is shown graphically in Figs. 24 and 25. 











A METHOD OF SOLUTION OF FIELD PROBLEMS BY MEANS 
OF OVERLAPPING REGIONS* 


BY 


H. PORITSKY anv M. H. BLEWETT 
General Electric Company 


1. Introduction. In problems involving the determination of fields, it often hap- 
pens that the region R for which the field is to be determined is difficult to handle 
directly, but can be broken up into several overlapping regions Ri, Re, - -- for each 
of which the field can be determined by standard methods. We suppose that the break- 
ing up is carried out in such a manner that every point of the region R falls into at 
least one of the regions Ri, Re, - - - . In the following, Schwartz’ “alternating proced- 
ure” is applied to the solution of field problems for such regions R. 

To illustrate, let us consider the determination of a function u harmonic over the 
region R shown in Fig. 1, bounded by two circular arcs ABC and CDA with centers 
at O and O’. For simplicity we assume that the radii of the two circles are equal and 
the boundary values of u are symmetric about the straight line through A and C. 
It will be noticed that by completing the circular arcs by means of the dotted curves 
AEC and CFA one obtains the circular regions over which the determination of a 
harmonic function in terms of boundary values is well known. Here R is the region 
bounded by the solid circular arcs ABC and CDA, 
while the regions R; and R; are the circular regions 
bounded by the complete circles with centers at O 
and O’. The problem then is to utilize the relatively 
easy solution of the Dirichlet problem for the circu- 
lar regions R; and R: in an efficient manner toward 
the solution of the Dirichlet problem over R. 

This is done by assuming the values of u over 
the arc A FC; the solution of the Dirichlet problem 
for the circle R; with center O, based on these assumed values and on the known 
boundary values over ABC, is then utilized to furnish the values of u over AEC. 
The procedure is then repeated by solving the Dirichlet problem for the circle R2 
with center at O’, and the values over A FC are recalculated. By alternating between 
Rand R’ in this way, continual improvement of the values of u over both arcs A FC 
and AEC results; in the limit this leads to a solution of the Dirichlet problem for the 
region of Fig. 1. 

In the following we shall illustrate the procedure, not for the Laplace equation 


Vu = 0, (1.1) 





Fie. 1. 


but for the equation 
(V? + k*)u = 0 (1.2) 


which is encountered in wave motion under the assumption of sinusoidal oscillations, 
for the region shown in Fig. 2. Other cases of interest in connection with (1.2) which 
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can also be treated by the present method are given by the “open end correction of an 
organ pipe,” wave passage through a change of width of a channel, 7-sections, etc. 

2. Wave propagation around a corner. We consider a solution of the differential 
equation (1.2) for the region shown in Fig. 2; this solution is to satisfy the boundary 


conditions 


Ou 
—=0 on DOG, EBF, (2.1) 
on 
u = A,e'*= + Bye-*** for large positive x, (2.2) 
u = Az3e**y + Bse-**y for large positive y, (2.3) 


where A;, B,, A; and B; are proper constants. Equations (2.2) and (2.3) can be de- 
scribed physically by the statement that u be- 
haves as a plane wave at infinity. 6 > F 

The above problem is encountered in the y| “ea 
propagation of a transverse electromagnetic . 
wave around a corner or through a channel the 
section of which is shown in Fig. 2. Here the © 








eee 








: 
channel has an infinite depth in the z-direction; ° los p | 
the field components are assumed to be inde- 4 : D 
pendent of z, and the only non-vanishing mag- ' | 
netic field component is H,. At the boundaries, } a - x o 
which are assumed to be metallic and perfectly 
conducting, the electric field is normal; this Fic. 2. 


leads to the vanishing of the normal derivative 
of H,, i.e., 0H,/8n=0. Formulation of the field in terms of H, leads to the above 
problem. 
On account of the vanishing of the normal derivative over the y-axis, reflection 
across it is possible, thus extending the region of 
6,4, Fig. 2 into the region shown in Fig. 3. This re- 
| t flection is carried out in accordance with the 
y relation 
r u(— x, y) = u(x, 9). (2.4) 


In view of this reflection the behavior of u at 
B' 8 x= -— © is given by the expression 


“= Bye*** 4. A,e7**, (2.5) 


ee ee 


= Bp" A, oP —> A, 





— Ao? B, <B, 
As a result of this reflection the semi-infinite 
strip DOCE of Fig. 2 can be replaced by the 


2-way infinite strip of Fig. 3 


° 
pe-- 
Pad 





Fic. 3. 
—-xox<2r< 2, 0<y<b. 


Similar reflection can be carried out across the lower boundary y=0 of Fig. 3; 
this allows us to replace the semi-infinite vertical channel by a vertical channel ex- 
tending to infinity both up and down. A proper behavior for u at y=— © can be 
obtained from (2.3). 
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The general procedure which was outlined in §1 is applied to the present case. First, 
we consider the strip 0<y<b of Fig. 3, and assume values for 0u/0n over the dotted 
part B’CB of its upper boundary. Since 0u/dn vanishes over the rest of its boundary 
and the behavior of u at © is described by (2.2) and (2.5), it is possible to determine 
u at any point interior to this strip. This determination is carried out by means of a 
Green’s function G. The derivation of G will be described presently; for the present it 
will suffice to say that the value of u at an interior point P of the strip is given by the 
relation 

1 > /du 
Up = u(xXo, Yo) = 2B, cos kx + — (<*) Gdx. (2.6) 
2 ~b \OY/ yao 
G has a pole at P=(o, yo), and (2.6) requires that G be evaluated on the dotted line 
B’'CB. After u is determined in this way, differentiation of (2.6) with respect to x 
allows one to compute 0u/dx, and in particular to determine this derivative over AB. 
Turning now to the infinite vertical strip 0<x <b, we repeat the same procedure and 
determine the function u at any point interior to this strip; in particular, we evaluate 
u and du/dy over CB. The process is then repeated. 

The definition of the Green’s function for the differential equation (1.2) and the 

boundary condition (2.1) for a finite region R is specified by the following: 


a) G satisfies (1.2) everywhere in R except at the pole P; 
b) 0G/dn vanishes along the boundary of R; (2 
c) at the pole P, G becomes infinite like —1n r’, where 7’ is the distance from P. } 


We apply Green’s theorem in the form 


Ov Ou 
f [u(V? + k2)v — (V2 + k)uldA = se — 9 ~) ds (2.8) 


On on 


to the region R, exclude the neighborhood of the point by means of a small circle of 
radius ¢ and let € approach zero. This yields the equation 
Be (2.9) 
“yp = — —Gds, i 
2r On 
where the integration is carried out over the boundary of R. In the present case, for 
the infinite strip 0<y<b special additional considerations are required. It will be 
assumed that in addition to the requirements (2.7) the Green’s function G behaves at 
infinity like a divergent plane wave. Solutions of (1.2) which depend on x only are 


et ike, (2.10) 
We consider the wave equation 
—— = ¢*V?U, (2.11) 
and look for solutions of the form ue**!. If we set k=w/c, we find that wu satisfies 


Eq. (1.2), and that e*** represents a plane wave traveling in the direction of positive x 
q I I 
while e~*** represents a similar wave traveling in the direction of negative x. It will 
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be assumed that at x= +o the Green’s function G behaves like a divergent plane 


wave of the same amplitude atx=+ © asatx=—o. 
It will be assumed that the dimension 6 satisfies the inequality 
b < a/k. (2.12) 


Physically this assumption means that the width b of the strip is less than half the 
wave length \/2 =7/k of a plane wave at the frequency in question. The effect of this 
assumption and the features which arise when it is not satisfied will appear presently. 

First, we place the pole P on the y-axis. We shall obtain G as a series in the form 


°c 
G = > Ants (2.13) 
n=0 
where u, are product solutions of the wave equation (1.2), i.e., u, consist of the prod- 
uct of a function of x and a function of y; more explicitly, 


uo = exp [ik- | x | |, 


ny 4 nn a4 " (2.14) 
Un = COS ; exp | (=) k =|], (n > 0). 


These product solutions u, (n>0) have been chosen so that they don’t become infinite 
at x= +, while uw» represents a divergent plane wave. If the inequality (2.12) were 
not satisfied, several radicals in u, (n>0) would be imaginary, infinitely large values 
of u, could not be avoided, and additional stipulations regarding the behavior of G 
at infinity would have to be made. 

The functions u, are symmetric about the vertical line x =0 through the pole P, 
and continuous at x =0. However, 0u,/0x is discontinuous at x =0, the discontinuity 
being 








r— ik for n= 0, 
Ou, \ 
Al - = < nx \? ny > (2.15) 
Ox 2 —)}-— k*cos for n> 0. ) 
k b 


Thus each term u, may be considered as representing the wave function correspond- 
ing to a sinusoidal distribution of sources* over the line x =0. The density o of the 
sources is given by the familiar condition from potential theory 








Ou 
discontinuity in normal derivative = a(=*) = — 2re, (2.16) 
x 


and in the present case is given by 


1 n*x? nTy 
oo, = — — k? cos ; . (2.17) 


ig 5b? 














* By a “source” is meant here a solution of (1.1) which depends only upon the distance r from a 
fixed point, is singular at r=0 like —In 7, and behaves at infinity like a divergent cylindrical wave.The dis- 
tributions of such sources satisfy continuity-discontinuity relations similar to those in the case of log- 
arithmic potentials. 
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For (2.13) this yields 





/n* 3? nny 


— k*? cos ——-: 2.18) 
b? ( 








: a. 
o=— DAs f/ 
T n=0 
Let us now consider the concentrated point source at the pole P, and express it as a 
Fourier series of cosines over the interval x =0, 0<y <b, obtaining 
1 r nT Yo nTy 
o = —+— )> cos — cos — (2.19) 
b P ant b b 
where x =0, y=yo are the coordinates of the pole P. Solving for A,, we obtain for 
the Green’s function G the Fourier series 





Ss) "Ac / ’ / 2 
+> = > : - COS — exp | / (=) — k?.| x| , (2.20) 
ant V(nx/b)? — BP b b f 

Due to the behavior of G at infinity it is found that after applying the Green’s 
theorem (2.8) over the rectangular region —/’<x</ and letting / and /’ recede to 
infinity, certain additional terms R’ and R” arise from the boundaries x =/ and x =I’. 
Equation (2.9) is now replaced by 





1 1 Ou(x, y) : 
u(xo, Yo) = —{ ——| Gdz+ R + R”, (2.21) 
2rd yp oy y=b 
where 
1 > /OUu 0G 
R’ = — (— G-—u )ay ; (2.22) 
2rJo \Ox Ox z=l 
i 7 Ou dG ; 
R” = — —-—G+u4— )ay : (2.23) 
2rd o Ox Ox ee 


In view of (2.4), (2.2), (2.20), (2.21) and (2.22), Eq. (2.21) can be given the form (2.6). 
As explained above, in the present case not (2.6) but its x-derivative will be found 
useful. Differentiation of (2.6) yields 


du(xo, Yo) : 1 > /du 0G 
es = — 2kB, sin kx + — a — dx. (2.24) 


OX 47 / _b oy y=b OXo 





To obtain this equation, the integral in (2.6) has been differentiated under the in- 
tegral sign; this is permissable since the limits of integration are independent of xo. 
Since (0u/dy)y-» is also independent of xo, only G has to be differentiated. The ex- 
plicit form of (2.24) is given by the relation 


Ou 1 jee? 1p 2 
g(Yo) = (=) = — 2RkB, sin kb — -f I x) | K - be K. | dx (2. 25) 
OXo Zo=b b b = 
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where 


beet Ou . 
f(x) = (~) : Ky = } exp [— ik(x — 6)], 
oy y=b 


‘ NT Yo ni \? 
K, = (— 1)" cos —— exp| 4/(“*) - B(z- s) | 
b b 
A similar expression holds for (0u/dyo) along CB; 


Ou : 1 b 
f(%0) = () = — 2kB; sin kb — —f g(y)[Ko+ >> Kaldy = (2.27) 
vob —b 


re) Yo 


(2. 26) 





where Ko, K,, are as in (2.26) but with the coordinates interchanged. 

In applying (2.25) one must assume not only f(x) but also B;. Likewise in applying 
(2.27), Bs is required along with g(y). Furthermore, A; and A; are essential to the 
complete solution. In this connection it is advisable to keep the following relations 
between f(x), g(y) and the constants A:, B;, As and B; in mind: 





1 b 
A, =B,- — { xje**zdx, 2.28 
mun A? _J i, 
1 b 
A; = B; —- — f y)et*udy, 2.29 
3 3 ikb 8 Y ( ) 
b 
Arett® — Bye-ik> = f y)dy, 2.30 
1 1 Dikb fo j ( ) 
1 b 
Asei*® — By = —— f x)dx. 2.31 
3 3 rT _J# ( ) 


These relations enable one to express A:, Bi, A3 and B; in terms of f(x) and g(y). 

The relation (2.28) is established by applying (2.6) to u(xo, yo) for xo so large 
that G reduces to its first term in (2 20), and comparing the result with (2.2). A simi- 
lar derivation over 0<x<b yields (2.29). As regards (2.30) it is established by ex- 
panding 0u/dx in the horizontal strip 0<y<b in a series of cosines of mmy/b and 
comparing for large positive x this expansion with 0u/0x as derived from (2.2); a 
similar procedure applied over the vertical strip 0<x<b to du/dy leads to (2.31). 

In the present example, in view of the geometric symmetry of the region of Fig. 2 
about the diagonal OB, any function u over the region can be expressed as the sum 
of a function which is odd about this diagonal, and one which is even about it. The 
calculations outlined are simplified considerably for even and odd functions u, are 
quite similar for the two cases and will be illustrated for the odd case. 


In the odd case, 


A3 =e Aj, B; eoo- B,, (2.32) 
g(x) = — f(x), (2.33) 
and the integral relations (2.28)—(2.31) reduce to 
1 b ro 
g(y) = — 2kB, sin kb — 7 f(x) | Ke+ >» K.| dx, (2.34) 
a ' n=] 
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1 b 
ee — | (x)e-#zdx, 
Ree 2ikb _ ) 
ea (2.35) 
A,e*® — Bye**> = — f (x)dx. 
; 2ikb _J ) 


After (2.25) has been applied for an initially assumed f(x) curve, the resulting g(y) 
shape, changed in sign and plotted against x, can be considered as the next approxima- 
tion to f(x), in view of the relation (2.33) and the symmetry of the region about x =y. 

From (2.35) the coefficients A; and B; may be determined in terms of 
(0u/dy)y-.=f(x). This yields 


1 b 
A — eS 3 —— [ x eik(z+b) _ 1 1x, 
4id sin B _I )| le 
: (2.36) 
B, = — ———— f(x) [eto — 1] dx. 
4kb sin kb J _, } 


The procedure used consisted in assuming f(x), calculating Ai and B, from (2.36), 
then applying (2.34) to calculate g(y), and using the shape of the latter with the sign 



































TABLE 1 
—— — — ee —— ee - > 
Ye x | ¥=K. bG/2e Yo x | UK bG/2x 
= 
1b 9b | —.4347 | .0613+.0627: 5b —.1b | —.0011 | 0920+ .4912% 
.7b | —.3025 | .1625+. 1841: —.3b | —.0003 | —.0308+ .4991i 
5b | —.1934 | .2111+.2938% — .5b | O | —.1545+.4775% 
.3b | —.1151 2034+ .3852i — ) 0 | —.2686+.4217i 
.1b | —.0682 1448 + .4524i — .9b 0 — .3607 + .3410i 
—.1b | —.0392 | .0545+.49124 $$ | |__| —_—_______—_ 
—.3b | —.0222 | —.0527+.4991i 7b .9b | —.1433 | .3527+.0627i 
—.5b6 | —.0126 | —.1671+.4775i 7b | 1255 | .5905+.1841: 
—.76 | —.0070 | —.2756+.4217i ‘sb | .1149|  .5194+.2938% 
—.95 | —.0040 | —.3567+.3410i 3b | .0725 | .39104+.3852i 
-- 1b | 0429 | .2559-+.4524: 
.3b .9b | —.4143 | .0817+.0627i ~, $5 0244 | .1181+.4912% 
.7b | —.2546 | .2104+.1841i a 0138 | —.0167+.4991 
5b | —.1458 | = .2587 + .2938: — .5b .0078 | —.1467+.4775% 
.3b | —.0806 | .2379-+.3852i —.7b .0044 | —.2642+ .4127i 
.1b | —.0453 | .1677+.4524i —.95 | .0025 | —.3582+.3410i 
—.1b | —.0252 0679+ .491 23 
—.3b | —.140 | —.0445+.49914 .9b .9b | 1.1279 | 1.6239+.0627: 
—.5b | —.0078 | —.1625+ .4775i .7b .5691 | 1.0341+.1841i 
—.7b | —.0044 | —.2730+.4217i 5b | .2862 | .6727+.2938i 
—.95 | —.0025 | —.3632+.3410i 3b | .1387 | .4572+.3852i 
—_——— 1b | .0744 |) .28744+.45247 
.5b .9b | —.3348 | 1612+ .0627 —.1b | .0410 1341+ .4912 
.76 | —.1373 | .3277+.1841i —.3b| .0226 | —.0079+4.4991i 
.5b | —.0437 | .3608 + . 2938 —.5b| .0126 | —.1419+.4775i 
.3b | —.0133 | .3052+.3852% —.7b | .0070 | —.2616+.4217% 


1b | —.0039 | .2091+4 .4524% —.9b | .0040 | —.3567+.3410i 
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changed as the starting point of the next step. To prevent the solution from becoming 
infinite, at each step f(x) is divided by A1, thus yielding the case A; =1. In the follow- 
ing numerical work the assumption )=0.2A, kb=72° is made. 

Although from physical considerations one would be able to make a reasonably 
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ai 5 & (20) 
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“1.0 -.8 “6 -.4 -.2 ° 2 4 6 rr) 1o- 
Fic. 4. 


good guess for the value of f(x), it was felt that in order to test the method thoroughly, 
the assumption 

Ou 

along BC, f(x) =| — = constant = 1, (2.37) 

DY] ymd 
would be more advisable. Making this assumption, solving for A; and B, from (2.36) 
(with b=0.2\, kb =72°), and dividing by Ai, we obtain 

bk — e*** sin kb 


A, = 1, B, = - 5 = .0634 — .999:i, 
bk — e~**> sin kb 





(2.38) 
f(x) 


Ou : 
~) = k(1.320 — 1.2382). 
OY) ymd 

The Green’s function was evaluated for five positive and five negative values of x, 
and for five values of yo, as shown in Table 1. The real part of this family of curves is 
shown in Fig. 4, the imaginary part being merely (2/6) sin k(x —b). These values, with 
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(2.38), were inserted in (2.34), the integration being made graphically with areas 
found by the trapezoidal rule, except near y=b. 

As y approaches 3), the value of g(y) increases so rapidly that extrapolation for the 
curve and the resulting graphical integration is difficult in this region. From physical 
considerations based on the fact that in a region which is small compared to a wave- 
length the function u behaves like a harmonic function, it may be shown that a fairly 
accurate approximation is obtained by assuming g(y) to vary as (y—b)~"? as y ap- 
proaches b. By picking two points y; and ye, two constants A and B can be found 
such that g(y) =A+B(b—y)~-" is fitted to the curve already drawn in this neighbor- 


hood for y <.9b; then the area is equal to 


b 3B 
J g(y)dy = A(b — y2) + — (b — yz)?/*. 

v2 2 
The resulting first approximation for 0u/dx is shown in Table 2. By means of (2.36) 
the values of A; and B; (the negatives of A, 
and B,) corresponding to these values were 
found to be Bz; =.1728 — 1.0062, A; =.998 —.1027; 











TABLE 2. 
— thus B,/A;=.2735 —.9801. 

g(y) = du/dy = (du/dx for In order to keep A; fixed at the value unity 
_ corresponding values of ~) which we have assumed, we retain this value of 
e —[1.004— .9445i] the Bi/A1 ratio and rename it B, as before. We 
3b —2[1.054— .9924] must then divide the values of du/dy in Table 2 
Sb —k[1.173—1.103:] by A;. Reinsertion now into (2.27) gives us the 
7b 6 | «= —k[1.382—1.299:] second approximation to 0u/dy shown in Table 
-9b —k[1.382—1.2994] 3. The corresponding A; and B, yield the ratio 
-% —k[1.876—1.7614] B,/A,=.2658 —.9607. The third approximation 





is then carried out in similar fashion, with the 
the results shown in Table 4. In this case, we 
have B,/A;=.266—.964i. The approximations to 0u/dy are shown in Fig. 5. Figure 6 
shows the ratio B,/A; and thus we see that this ratio is converging toward the value 
.266 — .9627, with the absolute value .997. 


TABLE 3. TABLE 4. 


The second approximation. The third approximation. 











f(x) =du/dy 








x f(x) =du/dy x 
.1b k[1.119— .804:] 1b k[1.081— .785:] 
.3b k[1.176— .844i] .3b k[1.132— .821i] 
.5b k[1.309— .960i] 5d k[1.249— .911:] 
.7b k[1.562—1.111:] .7b k[1.515—1.095:] 
.9b k[2.100—1.619:] .9b k[2.12 —1.557i] 








A similar calculation could be carried out for a function which is even about the 
diagonal OB. The results of this, together with those already found for the odd func- 
tion, would enable us to cover all cases involving a corner with these dimensions. 
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3. An alternative method of procedure. The procedure described and illustrated 
in the preceding sections can also be applied in a different way. Basically, the calcula- 
tion was carried out by first assuming the field over the line CB in Fig. 2, then cal- 
culating it over the line BA. It is possible to carry out the same calculation by assum- 
ing the field over CB not as a function of x or as a curve, but as a Fourier cosine series 


in x, 
nw x 





f(x) = ae >> Ca cos . (3.1) 


dy y=b 
Similarly, g(y) =@u/dx over AB can be converted into a similar Fourier cosine series 
in y, 
Ou | 
s(y)=—| = 2D, cos 


Ox z=b 


ny 


(3.2) 





Applying (2.34), (2.35) and (2.26) to the calculation of g(y) from f(x), we obtain 


MTYo 


g(yo) = >, Dn cos — —. 


b 


nn x 
exp [— ik(x — b)] cos dx 
b 





1 
= — 2kB, sin kb — — Mi 
,sin bY 


n ~ 


: ee... nTryo (° nNwx //m naa 
“—— > d (- 1)"C cos = cos ; exp / (= — (kk — i) as (3.3) 
m=1 


n=0 —-b 








This leads to integrals involving a cosine and an exponential in z. After these integra- 
tions are carried out, each one of the coefficients D, of the expansion (3.2) turns out 
to be linearly dependent upon the coefficients C,. Thus, instead of being given a curve 
f(x) and computing from it the curve g(y), one starts with B; and a series of coeffi- 
cients C, represented by the Fourier expansion (3.1) and ends up with the coefficients 
D, by applying (3.4). The explicit relation between these two sets of coefficients is 


Dy = — 2kB, sin kb ++ >) PorCn, Dm = >> PmiCn for m> 0, (3.4) 


n=0 
where 


(— 1)"ikb(1 — ei) 
Po, = ——— me 


2(n?x? — kb?) 
(— 1)"+™+1(m2x? — R%)?)'/2}1 — exp [— 2(m?x? — k%*)1/2]} 


Pu © ————— m > 0. 
(m? + n?)x? — kb? 











The matrix P,»,, thus takes the place of the series of curves 0G/0x which were given in 
Table 1 and shown in Fig. 4. A similar set of equations expresses C, in terms of D, 
and B3. 

By proceeding as in §2 with the field which is odd about the 45° diagonal OB, it is 
clear that for the final field the coefficients D, should be the negatives of the coeffi- 
cients C,. For the individual successive approximations, this of course is not neces- 
sarily the case. The calculation of the next improvement can be carried out by start- 
ing with D,, changing their signs and putting them in place of C, in (3.4). 
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It is possible to replace Dn by C,, in (3.4). A solution of the resulting equations 
would lead to a complete determination of the field problem. However, the solution 
of the resulting equations itself involves some method of successive approximation ; 
hence, this procedure is not advisable, and the successive calculation of C’s and D’s 
appears to be preferable, since it agrees in spirit with the method outlined above and 
constitutes just a variation of it. 











SOLUTION OF LINEAR AND SLIGHTLY NONLINEAR 
DIFFERENTIAL EQUATIONS* 


BY 
S. A. SCHELKUNOFF 
Bell Telephone Laboratories 


Considering the practical importance of linear differential equations of the second 
order, or the equivalent systems of the first order equations, it is surprising that trea- 
tises give little attention to effective and sufficiently general methods for their solu- 
tion. The treatises seem to be concerned primarily with power series expansions, 
Picard’s method of successive approximations, numerical methods based on difference 
equations—methods which in theory are applicable to almost any differential equa- 
tion and which are practically useless in the case of wave equations. On the positive 
side, in treatises on mathematical physics one finds a very effective asymptotic ap- 
proximation which in this country is known as the Wentzel-Kramers- Brillouin ap- 
proximation and in England as Jeffries’ approximation and, of course, the Rayleigh- 
Schrédinger perturbation method. The former has its obvious limitations and the 
latter is suitable only for a special class of boundary value problems. 

Our purpose is to call attention to another perturbation method which we de- 
veloped several years ago in connection with the antenna problem. As time went on 
the virtues of the method became increasingly apparent. Searching for previous 
references to this method, we came across one by Bécher! to a paper by Liouville.” 
In Liouville’s paper we have found the Jeffries-Wentzel-Kramers-Brillouin approxi- 
mation and a thorough discussion of the usual boundary value problem and associ- 
ated orthogonal series but very little that has any direct bearing on the present paper. 

The method is based on the idea that solutions of linear differential equations 
may be regarded as distorted or “perturbed” sinusoidal or exponential functions—the 
same idea which is back of the asymptotic approximation, of the Rayleigh-Schré- 
dinger method, and of the Sturmian theory. It is hardly surprising that this method 
gives better results than Picard’s method which regards the solutions as perturbed 
straight lines; but the difference is so remarkable that it deserves a special display in 
a separate note. In this paper, we restrict ourselves to an outline of the procedure 
and a statement of specific formulas reduced to a point where only simple integra- 
tions are needed in any special case. The exposition is based on the second order equa- 
tion; the extension to higher order linear equations is simple enough. When it comes 
to nonlinear equations, excepting those which are only slightly nonlinear,f the virtues 
of the method are not quite clear at present. There is no question that the results 


* Received July 6, 1945. 

1 Maxime Bécher, An introduction to the study of integral equations, Cambridge University Press, 
Cambridge, 1914. 

2 Joseph Liouville, Mémoires sur le développement des fonctions ou parties de fonctions en séries dont les 
divers termes sont assujettis a satisfaire d une méme équation différentielle du second ordre, contenant un 
paramétre variable, J. de Math., 2, 16-35, 418-436 (1837). 

t The meaning of “slightly” depends on the goodness of results expected from the process. Beyond 
that we shall not attempt to define it. 
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should be better when compared to those obtained by Picard’s method; but the more 
complicated technique for numerical calculations may offset the advantages. This is 
something to be explored. 

Suppose that our problem is to find the solutions of 


dV dI 
——— — Z(x)I, —_—_= — Y(x)V, (1) 
dx dx 
subject to the initial conditions 
V=V(a), T=T1(a), if x=a. (2) 
Picard simply integrates (1) and obtains a pair of integral equations 
Viz) = Va) — f z@rea, 12) = 1@) - f VOV@ude. (3) 


Thus the stage is set for successive approximations and the solution is obtained 


in the form of the infinite series 


V(x) = Vo(x) + Vilx) + Vole) +---, I(x) = Io(x) + I(x) + Ie(x) +--+, (4) 


where 
Vo(x) = V(a), V,,(2x) = -f Z()In-s(§) dE, 
¥ (5) 


To(x) = I(a), T,(x) = -f Y (&)Vn_a(€)dé. 


This procedure is so simple that it would be easy to overlook the fact that in sub- 
stance we are regarding the solutions of (1) as perturbations of the solutions of 


dV dI 

— = 0, —=0, (6) 

dx dx 
and that we are dealing with a special application of a much more general perturba- 


tion method. Let* 


Z(x) =Zo(x)+Z(x), Y(x) = Yo(x) + (a), (7) 
and suppose that the solutions of 
dV» dIo . 
—_— = — Zo(x)To, i Vo(x)Vo, (8) 
dx dx 


subject to the initial conditions (2), are known. Then the solutions of (1) are identical 
with those of the following integral equations 


V(x) = Volz) — f “Z(OUOVile, Hak — f P(e) V(E)V a(x, Bak, 
a a (9) 
I(x) — f Z(O1(OT1 (x, Edt — f V(é)V(8)I2 (x, &)dé, 


* In substance the theorem implied by equations (7), (8) and (9) is hardly new; but we have been 
unable to find its statement in just that form. 


I(x) 
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where V1(x, £), Ti(x, &); Ve(x, &), Ie(x, &) satisfy (8) and are subject to the following 
conditions 


Vig) =1, W1(&,§) = 0; Valt,t)=0, 2(, &) = 1. (10) 


Essentially the procedure is to regard —Z(x)I(x) and — Y(x) V(x) as known functions 
and to write the general solution of the corresponding nonhomogeneous linear equa- 
tion. The verification of the identity of the solutions of (1) and (9) is perfectly 
straightforward. If Ioi(x), Io2(x), are two linearly independent solutions of (8); then, 
as the reader can readily verify, ZoiJ02—JoeJ0, differs from Yo only by a constant 
factor. Bearing this in mind, we have 


T'o1(x)To2(&) — T’o2(x)Zox(€) 














Vi(x, £) = ———— nay 
T'o1(%)Io2(%) — I'o2(x)Io1(x) (11) 
To. To2 — Io2(x)Io 
SO a MR a 
I’ o1(x)To2(x) _ I’ o2(x)Io1(x) 
Similarly, 
Vo1(x) Vo — Voo(x)Vo1 
Vi(x, #) = —Zo(3) 1(x) Vo2(€) : (x) Vor(€) 
V'o1(x)Vo2(x) — V “o2( x) Vor(x) (12) 
T(x, §) = V'on(«) Vo2l(é) — V'o2(x)Vor(€) i 


Voi 2) Vox) — V'oa(x)Vor(x) 


Substituting Vo(x), Jo(x) in the integrands of (9), we obtain V;(x), [:(x); continu- 
ing the process we obtain solutions in the form (4). 

In Picard’s method Zo(x) = Yo(x) =0, which is the simplest possible choice. Natu- 
rally, the method will work well when Z(x) and Y(x) are small; otherwise it is far 
better to regard Zo(x) and Yo(x) merely as constants. If we are concerned with a finite 
interval, these constants may be chosen as some mean values* of Z(x) and Y(x)—the 
average values, for example; then for a=0 (9) become 


V(x) = Vo(x) -f Z(#)I(£) cosh To(x — £)dé + Kf Y(#)V(é) sinh T'o(x — )dé, 
0 (13) 
I(x) = I(x) + —{ Z(£)I(€) sinh T(x — dé -f Y(é)V(€) cosh I'o(x — é)dé 
0 0 0 
where 
Vo(x) = Vo cosh Tox — Kolo sinh x, To = VZ0l 0, Ko = V2Z0/Y 0, 
(14) 


Vo ; 

I)(x) = - x. sinh I'px + Jo cosh lox, Vo = V,(0), Io = I,(0). 
0 

In practice it is found that these equations represent a great improvement on 

Picard’s method and yet the integrations which have to be performed are not more 

difficult. If Z(x) and Y(x) are constants, Picard’s method leads to power series—not 


* Assuming that Z and Y do not change signs; if they do, it is best (although by no means necessary) 
to subdivide the interval. 
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a satisfactory form for wave functions. John R. Carson* employed Picard’s method 
for approximate solution when Z(x) and Y(x) are slowly varying functions and suc- 
ceeded in summing the series and obtaining the first order correction terms in a usable 
form; but any attempt to get the higher order terms by this method would seem to 
be out of the question. Theoretically, we should select Zo(x) and Yo(x) as near as 
possible to Z(x) and Y(x), subject to our ability to solve (8); but the integrations will 
be difficult to perform.* Thus we come back to (13) as the best compromise and it 
works very well. 
In the more explicit form the first order correction terms are 


Vi(x) = Vo[B(x) cosh Tox — A(x) sinh Tox + C(x) sinh Tox] 
— Kolo[A(x) cosh ox — B(x) sinh Tox + C(x) cosh Tox], 





Vo ; (15) 
I(x) = - = [B(x) sinh Tox — A(x) cosh Tox + C(x) cosh Tox] 
0 
+ Io{A(x) sinh Tox — B(x) cosh Tox + C(x) sinh Tox], 
where : 
1 (x) ~f lz KP | h 2M ofdé 
A(x) = — - cos ‘ 
x 7 ™ | Ky 0 fivaaé 
1 “rT Z . 
B(x) = =f ee Kt | sinh 21 oédé, (16) 
2 eo 0 
C(x) ~f [+ we? Jar 
xs) =— —— le. 
2JeL ke 





In some instances it is preferable to express the results in terms of progressive 
waves; then V(x) = Vo(x)+ Vi(x) and I(x) =Ip(x) +(x) become 


V+(x) = Kolo le Tes C(x)eTo _ E(x)eT*], 
I+(x) = Tit [ete — C(x)eTor + E(x)eh*]; 





(17) 
V-(x) = — Kola [eP* + C(x)efe? + D(x)eT], 
I-(x) = Ie [ePo= + C(x)ePo= — D(x]; 
where i 
1 “TZ a 
D(x) = A(x) + B(x) = —f | ~“ Ko? | e*Totde, 
2 0 Ko (18) 
E(x) = A(x) — B(x) f'[z Ket | —2,08ge 
WX) = ALL “= ) , Ko 0 é ‘rT ° 


Equations (14) and (15) express the solutions in terms of V and J at the beginning 
of a finite interval (0, /); one also often wants the corresponding expressions in terms 
of the final values. These are 


? John R. Carson, Propagation of periodic currents over nonuniform lines, Electrician, 86, 272-273 


(1921). 
* This objection would not apply in strictly numerical handling of equations. 
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Vo(x) V(D) cosh To(l — x) + Kol (2) sinh To(/ — x), 


Vil 
I(x) = ~*. sinh I'o(l — x) + I(1) cosh To(d — x), 


0 


Vi(x) = V(D{ [B(x) — BD] cosh Po(d + x) — [A(x) — A(D] sinh Tod + x) 
— [C(x) — C(D] sinh To(l — x)} 
+ Kol(l){[B(x) — B(D] sinh To(l + x) — [A(x) — A()] cosh To( + x) 
— [C(x) —C(D)] cosh To(l — x)}. (19) 
V()) , 
I(x) = x [A (x) — A(D)] cosh To(l + x) — [B(x) — B()] sinh Po(/ + x) 


— [C(x) — C()] cosh To(l — x)} 
+ I(l){[A(x) — A()] sinh To(l + x) — [B(x) — B(D)] cosh Po(l + x) 
— [C(x) — C()] sinh To(l — x)}. 
Suppose now that the interval is infinite and that Z(x) and Y(x) are slowly vary- 
ing functions. In this case, there exists the Liouville-Jeffries-Wentzel-Kramers- 


Brillouin approximation 


ll 


+ hey EGET op | +f ‘ rat], 
7 (20) 


I(x) = AW K(x%0)/K(x) exp l= f rat]. 


0 


V (x) 


where 
K(x) = VZ(x)/Y(x), I(x) = VZ(x)¥(z). (21) 


To the communication engineer these approximations seem natural even without 
formal analysis. He would reason as follows. If the “characteristic impedance” K (x) 
is independent of x, a progressive wave moving either to the left or to the right would 
suffer no reflection; it is only the sudden changes in the impedance that causes reflec- 
tions. Hence the voltage V(x) and current J(x) associated with the progressive waves 
will be proportional to exp F Lf, P (x)dx)]. If K(x) is a slowly varying function, we 
can ignore the reflections and in the first approximation consider the line as con- 
tinuously “matched” and thus acting as a transformer. This means that the voltage 
will vary directly and the current inversely as the square root of the characteristic 
impedance: hence, equations (20). 

There are several formal derivations;* but the one which appeals to us most be- 
cause it corresponds closely to the physical argument is also the one which permits 
further improvements in the approximation. Let us consider the “transfer parame- 
ter” 0 


. d@ 
@ = f r()dt, —=T(2), (22) 
dx 


as the new independent variable. Substituting in (1), we obtain 


4 John C. Slater and Nathaniel H. Frank, Introduction to theoretical physics, McGraw-Hill Book Co., 
Inc., New York, p. 148 (1933); John C. Slater, Microwave transmission, McGraw-Hill Book Co., Inc., 
New York, p. 73 (1942). 
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dV dl V 
—=—K(@)I, —=-—~. (23) 
d@ dO K(@) 
Eliminating first J and then V we have 
v"(0) ert v'(0) -V=0 I”’(@) + c= I'(0) -IT=0 (24) 
K(@) ere K(@) rr 


If K(@) is constant, we have simple progressive waves as anticipated; otherwise, we 
introduce new dependent variables in conformity with our idea of voltage and current 
transformation 

V=[K@)]}}?V, I= [K@)}7. (25) 
Incidentally, this is the transformation which should remove the first derivatives from 
(24). Substituting, we obtain 

3(K’)? eV (K’)? yr 


v"@) = E +——_-—|v, 1’@)= E -— +> 


(26) 
4K? 2K 4K? 2K 





We now have not only equations (20) but also the quantitative criterion of their 
goodness: (K’/K)* and K’’/2K should be small compared with unity. 

To improve on (20), we could repeat the process beginning with (22); but the 
analytical work is simpler if we turn to equations (13) and apply them to an infinite 
interval, assuming of course that in the entire interval the bracketed quantities in 
equations (26) differ but little from unity. Thus, the solutions of 


d*y 





= [1 + f(x) ]y (27) 


9 
x? 


are also the solutions of 
v(2) = sol) + f @y@® sinh (x — Hat, (28) 
provided the integral is convergent. The solutions asymptotic to e** are 


y(x) > eF* + [toc sinh (x — é)dé, (29) 


or 


yt(a) ~e-® — fe f (8d + he f ee f(E) dé, 


ya)wetie f sodt—aer fo eyfeoae (30) 
From these equations we can obtain the well-known asymptotic expansions of Bessel 


functions as well as expansions of other types. 
The case in which @ =i8x, where B is a constant, occurs so frequently that a 


repetition is justified. Equations (26) become 


vay =-e74+(-2- |v, 7%) = eT +|— |r ep 
(x) = -—8 +| OK |.’ _ 2K 4K? 
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and the corresponding integral equations are 





(a) = Vala) +— fi Foam V(®) sin B(x — Oat 
x 0 4[K 2 2K(é) ei 
Te) =Tua) + =f [—* sO le sin B(x — dé 7 
x) = Io(x nie - _ - : é rp — E)aé. 
BJ.L2K(é 4[K(é]? 
Suppose, for example, that K(x) =Ko+kx; then, asymptotically, 
3ik : 
Vix) = + AV Ko t+ [1 —-— — | ere, 
88(Ky + kx) 
A tk 
I(x) = ———_ E + ——_——_ Je Fibs, (33) 
/Ko + kx 88(K ro + kx) 


In this case, however, the integrals in (32) can be evaluated in terms of sine and 
cosine integrals. Moreover, the complete result corresponds closely to the physical 
picture of reflection which invariably takes place when waves are traveling in trans- 
mission lines or media with variable characteristic impedance K(x). Thus 


V(x) = Al /Ko + kx e~#* + Ry\/Ko + kx e#*], 


e~ ib eisz (33’) 
sot wig eto) 


I(x) 
V Ko + kx V Ko + kx 


where Ry and R; are the first order reflection coefficients given by 


Ry = — 3Rr = — (3/4) exp (2iBk-'Ko) ci (2Bx + 2Bk-'Ko) 


~ 7 Si (28x + 2BkK,) + =| 


The succeeding correction terms represent successive reflections. The entire series re- 
sembles an asymptotic solution of the differential equation in question but it appears 
to be rapidly convergent. 

An another example, take the case of principal waves on a thin cylindrical antenna 


when 
. | 120 120 
K(x) = 120 log (2x/a), K'(x) =— K"(z) = ——- (34) 
x x? 
In this case we obtain 
V(2) = AVE 51 - Te | 3 
(x) = AV K(x - ———_—_——+} |e«~*=, 
26x 4flog ( (2x) var 2 log 2(x/a) (35) 
35 





ae | 
( 
I(x) A [1+ tf 4 1 I] . 
(xy) = — —— —- ee ee e7 Bz, 
/ K(x) 2Bx (4 flog sar 2 log (2x/a) 
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As the third example we shall take Rayleigh’s equation for a nonlinear oscillator® 
q+ (Rig + Rag*) + wg = 0. (36) 
By (13) we have 


1 t 
g(t) = go(t) — —f [Rig(r) a R3q*| sin w(t — r)dr, (37) 


where go(¢) is a sinusoidal function. If g=0 up to ¢=0, then go(t) =A sin wt. Substitut- 
ing in (37) and integrating, we obtain 


1 
g(t) = A sin wt — 4(Ry + 3w?R3A%)? sin wi— rte 3(cos wt — cos 3wt). (38) 
For a periodic solution we must have 
Ri + 3w*R3A* = 0; (39) 
then 1 
g(t) = A sin wt — op net 3(cos wt — cos 3wé). (40) 


Equation (39) is precisely Rayleigh’s equation for the amplitude of oscillations; equa- 
tion (40) differs from his equation in that ours contains a term proportional to cos wt. 
Our approximation satisfies the initial condition g(0)=0 while Rayleigh’s does not. 

Originally this work was undertaken to obtain convenient analytic approxima- 
tions to a number of problems in wave theory. It has since become apparent, however, 
that at least for a certain class of differential equations, the method would be suitable 
for numerical solution. The practicability of Picard’s method for this purpose has 
already been explored by Thornton C. Fry; the present method should be quicker. 
The rapidity of convergence will be discussed in a separate paper. 

5 Ph. LeCorbeiller, The nonlinear theory of the maintenance of oscillations, 1.E.E. Journal, 79, 361- 
378 (1936). 


6 Thornton C. Fry, The use of the integraph in the practical solution of differential equations by Picard’s 
method of successive approximations, Proc. 2d Internat. Cong. Math. Toronto, 2, 405-428 (1924). 
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A CYLINDER COOLING PROBLEM* 


BY 


SAMUEL A. SCHAAF 
University of California, Berkeley 


1. Introduction. The linear cooling problem for non-homogeneous solids has been 
investigated extensively by Rust,’ Churchill,? Carslaw,? Mersman,‘ and others. It is 
the purpose of this paper to obtain a solution for the corresponding cylindrical prob- 
lem. The method used is that of the Laplace Transform. 

2. The Problem. Let us consider an infinitely long circular cylinder of radius a 
and initial temperature To, instantaneously immersed in an infinite medium initially 
at zero temperature. Let the heat conductivities and diffusivities of the cylinder and 
external medium be respectively K, and h? (v=1, 2). Then if r is the distance from the 
axis of the cylinder and ¢ is the time, the following differential system is satisfied® by 


the temperature functions 7,/(r, t): 





2 fT) 1 oT; OT, / 
hy — —S} <= ose <a, >, (1) 
| Or? r or ot 
2 {°T 2 1 OT») OT» 
ho S$ ——_ + — - f = — r>e, > 6, (2) 
l ar? r or ot 
lim 7; = lm 7; > ©. (3) 
ta— rat 
OT, OT > 
lim K,—= lim Ke— i> 0, (4) 
ra— or rat or 
lim 7, = JT» Osr<a, (5) 
10 
lim 7, = 0 er (6) 
t—-0 


3. Solution. Let the Laplace transform of T,(r, t) be T,*(r, s), i.e., 
T(r, s) -{ eT (r, idt ce: 
0 


Applying this transform to (1)—(6), we obtain the corresponding set of ordinary dif- 
ferential equations containing s as a parameter; 


* Received June 18, 1945. 

1W. M. Rust, Jr., Integral equations and the cooling problem, Amer. J. Math. 54, 190-212 (1932). 

2. R. V. Churchill, A heat conduction problem, Philos. Mag. (7), 31, 81-87 (1941). 

3H. S. Carslaw, A simple application of the Laplace transformation, Philos. Mag. (7), 30, 414-417 


(1940). 
4W. A. Mersman, Heat conduction in an infinite composite solid, Bull. Amer. Math. Soc. 47, 956-964 


(1941). 
5H.S. Carslaw, Theory of heat conduction, Macmillan New York, ed. 2, 1921, Chapter I. 
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git gf BA T+sT# OSr< (1*) 
hy <§ ——_ + — —— > = — s <r<a, 
‘\ dr? r dr f : ; 7" 
2 (d?*T$# 1 dT# 
hy -4+ — } = sT} ¥> @, (2*) 
dr? r ar 
lim T* = lim TF; (3*) 
r+a— rat 
E adT#* dT# 
lim K,——= lim K,——- (4*) 
ra— ar rat dr 
The solution® of this system is 
To Ky /s)-I /s/h 
T#(r, s) = aE et See | (7) 
s D(vV/s) 
, J T 
T#(r, s) = To asl ¢ (ary s)-Kolrv s/he) (8) 
Ss D(vVs) 
where 
a, = a/h, a, = a/ho, a3 = Kyh2/Keh, (9) 
D(x) = asl) (a,x) K (ax) ated To(ayx)K¢ (ax). (10) 


The functions 7,(r, 4) may now be obtained by use of the complex inversion 


formula? 
ct+ik 


1 
T(r, 4) = lim — e*'T*(r, z)dz, c> 0: (11) 


Ave LATS cin 


In order to reduce these contour integrals to real integrals we must first establish 


the following lemma. 
LEMMA. D(z) does not vanish for | arg z| <r. 
Proor. We choose two numbers A and B, arbitrary except that 





0<A<K1<B. (12) 
Then since a, is positive (v=1, 2, 3), it will be sufficient to show that, when | arg 2| 
<4, D(z) does not vanish for any values of a, such that A Sa, $B (v=1, 2, 3). The 
proof now follows in four parts. 


i) There is a number R;>0, such that D(z) does not vanish for 
Asa, B, | arg z| < 4n, || < Ri. 


This is true because we may use the ordinary series expansions of the Bessel function® 


to write 


1 





Zz Qo 
D(z) =- 4 (ae = aa) > og + B(z), 


AF 
6G. N. Watson, Theory of Bessel functions, Cambridge University Press, Cambridge, ed. 2, 1944, 
p. 79. 
7D. V. Widder, The Laplace transformation, Princeton University Press, Princeton; Oxford Uni- 
versity Press, London H. Milford, 1941, p. 66. 
8 G. N. Watson, loc. cit., pp. 77, 80. 
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where B(z) is bounded in the neighborhood of z=0. The result is then evident. 
ai) There is a number R:>0, such that D(z) does not vanish for 


Asa 3B, | arg z| S 4m, |z| > Ro. 


To see this we use the well-known asymptotic formulae for these Bessel functions® 
to write 


e(al—a2)z 


a { (as + 1 + €1) _ (ag — ] a €)e~2a18} 


22V/ ajae 


D(z) = 


where €,—0 as z— © for k=1, 2, uniformly in the a, providing A Sa, $B (v=1, 2, 3). 
Clearly, D(z) can vanish only if 
az3—-1l+es. 


era — — 


- 6. 
agtit+ea 


But for sufficiently large | s|, say | z| >R:2, the right member is less than unity in 
absolute value. Hence for |z| >Re, this relation cannot hold with |arg z| <3. 
411) D(z) does not vanish for | arg z| =}. To see this, we let z=e!'*y (y real). 














Then!” 
Io(ed**y) = Jo(y), (13) 
Ko(elry) = — dix[Jo(y) — i¥o(y)]. (14) 
Hence 


D(yetir) = {asJé (ary) Vo(a2y) — Jo(ary)¥é (azy)} 
+ ifasJé (ai1y)J o(aey) = Jo(aiy)J¢ (ay) }. 
Therefore D(ye**) can vanish only if 
as] 9 (ary) Yo(ary) — Jo(ary) Vo (azy) = asJd (ary)Jo(ary) — Jo(aay)Jo (azy) = 0. 


But this is impossible since it would imply either the existence of a common root for 
at least two of these Bessel functions," or the vanishing of the Wronskian 





W [Jo(a2y), Yo(ary)] = 


Ta2y 


iv) We consider now the integral (see Fig. 1) 


ye 1 f D’(z) P 
ewer eels aa 


From 1), 7), and iz) it follows that D(z) does not vanish on C for all a, such that 
A Sa,=B (v=1, 2, 3). Now these Bessel functions are all analytic except possibly 
at z=0. Hence f(a, a2, a3) is continuous. 


* G. N. Watson, loc. cit., pp. 202, 203. 

10 G, N. Watson, loc. cit., pp. 77, 78. 

11 Tt is a well-known result that these Bessel functions have no common roots. See G. N. Watson, 
loc. cit., pp. 479, 480, 481. 
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Since D(z) has no singularities inside C, f(a1, a2, a3) gives the number of zeros of 
D(z) inside C. It can therefore take on only integral values; but this implies that 
f(a, @2, a3) is constant. 
Finally 
1 W’ [To(2), Ko(2) | 


1,1,1) = — z= 
JUL D = TJ Wits), Kels)] 





Hence f(a, a2, a3) =0, for all a, satisfying the relation A Sa,<B. Therefore D(z) has 
no roots inside C. Since the radii Rs and R, (see Fig. 1) are arbitrary, except that 























(¢ +i) 

6 } 

r\R3 Ra L, lan 

®, SS Ne 

l 
i 
|} (¢-id) 
Fic. 1. The contour C, consisting of the circu- Fic. 2. The contours /, T;, Ms, Ms, Z; and Ls. 
lar arcs |z| =R; and || =R, and the line seg- The radius of I; is p. 


ments on the imaginary axis joining them. The 
only restriction is that R;<R, and Ry>R2. 


R;<R, and Ry>Rz2, it follows that D(z) has no zeros in the entire right-half plane, 
which concludes the proof of the lemma. 

We now transform the contour integrals of (11) into real integrals. Let us consider 
Ti(r, t) first. According to the lemma just established, D(./z) does not vanish 
for larg V/z| S4r, i.e., for larg z| <7. Hence the integrand in (11) is analytic 
for |arg z| Sa, and we may (see Fig. 2) replace the integral along ] by the sum of 
the integrals over T,, Ts, T's, Zi, and Le. Using the asymptotic developments, we 


easily see that for large z, 


K6 (aav/z)To(rv/z /h1) 1 r x 
acne = 0 = exp} (— — a jVv zie. 
V2 hy 


D(V2) 





Therefore as \> ©, the integrals over , and I vanish, since ¢>0. 
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Near the origin, the term Ké (a2+/z)Jo(a1+/z) dominates the denominator, and 


hence 


. Ko (2/2) Io(r/2/ hi) _ 


iim ak 
ies D(vV/2) 





Therefore the integral over I’; vanishes with p (see Fig. 2). 
On Li, we set z=0°e**, ¢>0. Then, using (13) and (14), we obtain 


1 T, (8 ent 
rae e'T #(r, )ds=— f 
2rid 1, Tid vp o 
ae 


{i+ J o(ra/ ‘hy)V¢ 0 (a0) + Jo(ro/ha)Jo | (a2c) 
[asJ¢ (a10) V o(a2e) —Jo(ara) Vg (aoe) |+ilasJé law) J dest) — J o(ayo)J¢ (as) ] | 


On Lz, we set 2=07%e—**, ¢ >0, and obtain the conjugate of (15). Adding these and 
taking the limit as \> » and p— 0, we obtain finally, 

















~ sv © et Ti(ra/hy) Jo 
T(r, #) = a fae. (16) 
Ta 0 o? A(c) 
where 
A(c) = lasJ¢ (a10)Vo(a2a) — Jo(aic) Vo (2a) |? (17) 


+ [asJé (ay0)J (age) — Jo(ayo)J¢ (carga) ]?. 


Similarly 
Tas (7% e-* * Jd (are) [Jo(aic)C (a2, ra/h2) —as3J¢ (axc)C(age, ro/ he) | do, (18) 
— do, 





arene! Mane A(o) 

where 
C(x, y) = Jo(x)Yo(y) — Yo(x)Jo(y), (19) 
Cx(x, y) = Jo (x) Vo(y) — Yo (x)Jo(y). (20) 


These formulae constitute the solution of the differential system (1)—(6). 

4. Remarks. Since the Laplace Transform method is essentially a formal one, any 
solution obtained in this manner must always be verified. In the present case this is 
easily done.” 

It may also be shown, under certain conditions as to boundedness and continuity 
necessarily satisfied by any physical temperature distribution, that expressions (16) 
and (18) constitute the unique solution of the system (1)—(6).¥ 

In conclusion, the author would like to express his thanks to Professor G. C. Evans 
for his help in the preparation of this paper. 


2 For an example of the method see H. S. Carslaw and J. C. Jaeger, A problem in conduction of heat, 


Proc. Cambridge Philos. Soc. 35, 394-404 (1939). 


18 For an example of the method see W. M. Rust, Jr., loc. cit., p. 196. 
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ANALYSIS OF NUMERICAL SOLUTIONS OF TRANSIENT 
HEAT-FLOW PROBLEMS* 


BY 


CLARENCE M. FOWLER 
U. S. Naval Academy 


1. Introduction. The purpose of this paper is to present formal methods for es- 
tablishing the convergence of numerical solutions of transient heat-flow problems, 
and to derive expressions for these solutions in terms of the initial temperatures and 
boundary values. 

In general, heat-flow problems are classified under two groups, steady-state flow 
and transient flow. Steady-state problems are solved numerically by the relaxation 
method. Many papers dealing with the actual numerical work have been written, 
and Temple! has established the validity of the relaxation method under various 
boundary conditions. Moskovitz? has derived an expression in terms of the boundary 
temperatures for the steady-state numerical solution of a rectangular bar. 

Although considerable work has been done on the actual application of numerical 
methods to transient heat-flow problems, very little has been written about the’ prob- 
lems of convergence and the expression of solutions in terms of initial and boundary 
values.’ These last two considerations are the objects of this paper. 

Two restrictions which simplify the analysis are placed on the examples consid- 
ered here. First, only the one-dimensional slab is considered; secondly, the initial 
temperature distribution is assumed to be constant over the slab. However, by ex- 
tensions of the methods used, solutions of problems concerning two- and three-dimen- 
sional rectangular objects with arbitrary initial temperature distributions are readily 
derived. 

The various boundary conditions which have been studied include the following : 
the temperature at the boundary is given, and is either a constant or a function of 
time; the boundary is insulated; there is a constant energy input at the boundary; 
there is convection at the boundary. The author has made no attempt to consider all 
possible combinations of boundary conditions, but has tried to include enough repre- 
sentative cases to illustrate the methods. 

The procedure followed throughout the paper has been to consider each example 
as a whole, and to derive solutions of the problem and take up a study of the conver- 
gence, before proceeding to the next example. In some cases solutions have been ex- 
pressed in terms of a set of polynomials which are associated fundamentally with the 
difference equation; in other cases they have been expressed in finite Fourier series. 


* Received April 11, 1945. 

1G. Temple, The general theory of relaxation methods applied to linear systems, Proc. Roy. Soc. London 
(A), 169, 476-500 (1939). 

2 D. Moskovitz, The numerical solutions of Laplace's and Poisson's equations, Quart. Appl. Math. 2, 


148-163 (1944). 
* R. Courant, K. Friedrichs and H. Lewy, Uber die partiellen Differenzgleichungen der mathematischen 


Physik, Math. Ann. 100, 22-74 (1928). 
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Properties of the polynomials needed throughout the paper have been demonstrated 


in an appendix. 
2. The difference equation and boundary expressions. The basic one-dimensional 


difference equation satisfied by the numerical solutions is 





T o-1,1-1 + OT o.s~1 +T, 1,t—1 
Ti1= Saas (2.1) 
a+2 
where 
(Ax)? 
At = ——_—_—__ - (2.2) 
(a + 2)a 


Here, At is the time interval between successive time values, Ax is the distance be- 
tween successive points across the slab, ¢ is the time (in units of At), x is a space co- 
ordinate running across the slab (in units of Ax), T,,, is the temperature in the slab 
at time ¢ and position x, a is the thermal diffusivity of the material, and a is the 
modulus of the equation. 

As usually encountered, the difference equation has a=0 (Schmidt’s equation). 
Dusinberre* generalized Schmidt’s equation by introducing the modulus. The value 
of using an arbitrary modulus lies in being able to select an arbitrary time increment 
as well as space increment. This is not possible in Schmidt’s equation, since fixing Ax 
determines At. 

In dealing with convection, insulation, etc., at a boundary, it is always necessary 
to make some assumption to determine the numerical 
boundary expression. It should be emphasized that, for 
this reason there are several different expressions in use 
approximately representing the same boundary condi- 
tion. However, it is possible to consider only one of them 
here, which is deduced as follows. Figure 1 shows the 

Fic. 1. boundary (x=0) and the first two interior points of the 

slab. The boundary expression is derived by making a 

heat balance over the shaded half-segment. The heat gain by conduction throughout 
the time increment Af referred to the initial time instant t—1 is 


(Ti1-1 33) T 0,+-1) 


— hAAUT 94-1 — Ta) + RAAt ———— 
Ax 











? 


where / is the surface heat transfer coefficient, k is the thermal conductivity, 7. 
is the ambient temperature, and A is a unit area perpendicular to the slab cross sec- 
tion. This quantity is equated to the heat capacity gain 3(Ax)cpA(To,1:—To,:-1), 
where c and p are the specific heat and density of the material, respectively. For rapid 
surface cooling, it is necessary to keep Ax small, since it has been assumed above that 
To, will represent the temperature of the shaded segment. 
After equating the two heat quantities above, and simplifying, we obtain the 

boundary condition 

“ 2T 14-1 + (2 — 2N)To,1-1 + 2NT»o 

T 0,1 = — . ’ (2.3) 

a+2 








4G. M. Dusinberre, Numerical methods for transient heat flow, Trans. A.S.M.E. 67, 703-709 (1945). 
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where JN, the equivalent numerical transfer coefficient, is given by the relation 
N = hAx/k. (2.4) 
In the case where a boundary has a constant energy input per unit area, g, an 
analysis similar to that given above yields 
2714-1 + aT o,t-1 + 20 
To = ’ (2.5) 
a+2 


where 
QO = gAx/k. (2.6) 





For an insulated boundary, N=0 and (2.3) becomes 
272-1 + OT 01-1 
T 0,1 paca , 
a+2 


For simple boundary conditions such as temperature at x =0 held at uo, or tem- 
perature at x =/ fixed as a function of time f(t), the boundary conditions are simply 
T0,,=Uo or T;,,=f (2). 

3. Convergence. There are two distinct types of convergence to be considered 
here. The first type deals with the convergence of numerical solutions as the time 
becomes large. The second type shows that as the time and space increments At and 
Ax are allowed to approach zero, the numerical solutions become identical with the 





(2.7) 


corresponding analytic solutions. 

From an inspection of (2.1) and (2.2), it is seen that in applying numerical solu- 
tions to any particular example, there are apparently two arbitrary quantities, the 
space increment Ax, and the modulus a, which in turn determines the time incre- 
ment At. However, it is found from experience that if the value of a is taken too 
small, the calculated numerical answers oscillate and ultimately diverge as the time 
becomes large. The first type of convergence is concerned with developing criteria 
which impose a lower limit on allowable values of a which will then insure numerical 
convergence. Each example considered has such a criterion developed, since such 
criteria usually depend on the particular boundary conditions. Another related prob- 
lem is that of determining the steady-state distribution given numerically. It is shown 
that numerical solutions converge to the same steady-state values as those deter- 
mined analytically for the boundary conditions under consideration. 

Both of the problems discussed above pertain to actual numerical solutions where 
the space and time increments are finite, non-zero quantities. The second type of con- 
vergence is treated in §11, apart from the main body of the paper, since it does not 
deal with numerical solutions as applied, but rather to the limiting case where Ax 
and At approach zero. Under these conditions, the numerical solutions become identi- 
cal with the analytic solutions for all values of time and position throughout the slab. 

4. Particular solutions and contour integration. By substituting T,,,= F(t) sin 2x 
into (2.1), F(t) is found to satisfy the subsidiary difference equation F(é) 
=[(a+2 cos z)/(a+2)|F(t—1) and therefore F(t) =[(a+2 cos z)/(a+2)]*, from 
which we find as a particular solution to (2.1), 


a+ 2 cos 2\$ 
T24 = Gee sin 2x. (4.1) 
a+2 
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A similar analysis shows that sin zx may be replaced by cos 2x or e‘#*, 

Using (4.1) and extending it to the two or three-dimensional form, we can write 
down immediately solutions in Fourier series and integrals for rectangular objects, 
with arbitrary initial temperature distributions. However, such solutions are of little 
use as they converge too slowly. Instead of following the standard Fourier develop- 
ment, the author has found it expedient to consider only cases where the initial tem- 
perature over the slab is zero (which by a change in temperature origin includes any 
constant temperature distribution). This restriction, for the one-dimensional slab, 
allows the use of a method of contour integration which may be summarized as fol- 


lows. 
a) Having found that 


a+ 2cos2\'TA(z Biz 
Tea = (“= )| «) cos 2% + > sin 2x | 
a+2z Zz Z 


is a particular solution as long as A(z) and B(z) are independent of x and #, we 
formally integrate this solution with respect to z over the prescribed path (Fig. 2) 
in the complex plane. The solution is (4.2) below, where the functions A(z) and B(z) 
are determined so that (4.2) satisfies the boundary conditions of the problem. 











o 


1 a + 2 cos 2\ ee dz 
T2441 = —f (“——*") (4 (z) cos zx + B(z) sin zx )— - (4.2) 
Tid p a+2 z 


The path P is chosen parallel to the real axis, extending from + © to — ©, and is 
located a finite distance m above the real 














i axis, m being determined so that all 

fi 2 ie poles of the integrand lie below P. 
he < \ b) The integrand of (4.2) is shown to 
/ 2 % vanish over the arc, path R of Fig. 2, 
/ . E i \ except possibly at the slab boundary 
. Re r } points x =0 or x=/, when ¢ is given the 
ul ! | value zero. Then, as there are no poles 
1 ; “_ —" suanite ws enclosed by paths P and R, it follows 
‘ Q / from Cauchy’s theorem that at time 


zero T,,9=0, except possibly at the slab 
Fi. 2. boundaries. It follows that (4.2) is the 

solution to the problem, for it satisfies 

the difference equation, the initial condition T,,.=0, and the boundary conditions. 

c) The remaining step is the evaluation of the contour integral. This is accom- 
plished in one of two ways. In either case, the integrand of (4.2) is shown to vanish 
over the paths M and WN (Fig. 2) or over half these paths. 

1) For semi-infinite slab problems, the integral over P is evaluated in terms of an 
integral along the real axis and the residues of any poles lying between path P and 
the real axis. 

2) For finite slab problems, the functions A(z) and B(z) are generally such that 
the integrand is even-valued, and therefore the solution (4.2) may equally well be 
integrated over the path Q (Fig. 2) which is opposite path P. Thus, integration around 
the loop consisting of paths P and Q, and the paths M and WN shows that the required 
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integral is equal to half the sum of the residues at the poles enclosed by paths P, Q, 
M and N, since the paths M and N contribute nothing to the integral over the loop. 
All analysis of a purely rigorous nature has been omitted from the paper, but all 
doubtful cases have been tested for proper convergence and the vanishing of the in- 
tegrals over the paths outlined above. 
5. Semi-infinite slab. Boundary x=0 held at constant temperature 1, initial 
temperature zero. Let us consider the following equation 


u a+ 2 cos 2\! dz 
Ti = on (=o<"*) et? —. (5.1) 
a+2 


ni z 
To prove that this expression is the solution of the problem it is necessary to show 
that T,,9=0 and that 7o,.=%o. 

To show that 7,,9=0, we set t=0 and make the substitution z= Re‘*. We then 
integrate (uo/m) exp [ixR exp (ip)|dp over the path R, @ varying from 7 to 0. It fol- 
lows that as R- © this integral vanishes, except at x =0. Therefore, since there are 
no poles enclosed by paths P and R, it follows from Cauchy’s theorem that T,,o=0. 

To show that the integrand vanishes over the paths M and N, let us substitute 

= +R++iy where y varies from —m to +m, and let R-~ ; then it follows that these 
integrals vanish. To,;, the temperature at x =0, can equally well be integrated over 
path Q, since the resulting integrand is even-valued. Therefore, by Cauchy’s theorem, 
since there are no contributions from the paths M and N, To,1=%o. It follows that 
(5.1) is the solution to the problem, for it satisfies the initial and boundary conditions. 

To evaluate (5.1) it is convenient to integrate around the loop consisting of P, 
the half-path M, the real axis indented at the origin by a small semi-circle of radius e, 
and the half-path NV. In the work that follows, the real axis will be denoted by w to 
avoid confusion with the slab position variable. 

Since the integrals over M and N vanish, and there are no poles enclosed in the 


loop, 


uy (t*f/at+2cosw\! dw 
‘ —_ — CS) e7 zu —— 
wt a+2 w 
uo (f° fat 2 cos ce*\' @ 
+ ——_—_———— } exp lixe exp (ip) |do 
Wd a+2 
uo (t°/at+2cosw\' dw 
_— — Cao") et izw —___ = 
wt a+2 w 
Regrouping terms and letting e—0, we have as the solution 
2 £°%/at+2coswv\' . dw 
T 24 = ul = =f —————— } sin wx —— ]. (5.2) 
Tv 0 a a 2 Ww 


To express this solution in terms of the polynomials, we expand (a+2 cos w)' ac- 
cording to (12.7) and substitute into (5.2), 


Ta? uw] 1— par rs — — J [P.(t) + 2Pis(é) cos w + - 


dw 
+ 2Po(t) cos tw] sin wx =|. 
w 
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Combining the sine and cosine terms and making use of the identity (12.6a), 
P,_,(t) =P :4,(4), we have 





2 co r=+t : dw 
Tz, _ us| 1 _ (a + 2)* ‘ > Pri t(t) sin w(x + r) =] 
or 
1 = 5 
en = j 1— — }x Pritlt ’ Le 
melt oem Eat teat 6 


where the symbols {x+r} =1, 0, or —1 depending on whether x++r is greater 

than, equal to or less than zero. This notation arises from the fact that 

(2/2) f- (sin w(x+r)/w)dw= {x+r} in accordance with the above convention. 
From (5.2) it is easily shown that the solution converges as t—« when 


a2 0. (5.4) 


To show that the numerical solution converges to the analytic steady-state solu- 
tion for the same boundary conditions when t—>~, the following device is used: 
T.,, may be equated to the integrand over the real axis, and added to the sum of all 
residues of poles enclosed by the real axis and the path P. The real axis is to be in- 
dented with small semi-circles at all poles lying on it. In the steady-state value of 
the solution, the only contributions which remain as tf, are those which occur 
where z=2nz. All other contributions, including those along the real axis, drop out 
due to the rapidity with which the factor [(a+2 cos z)/(a+2) ]'-0 when z#2nz, as 
t— 

In the problem under consideration, there are no poles, and if we take into account 
the indentation at the origin, we find that the steady-state solution approaches wo. 
This is also the value given by the analytic solution of the same problem. 

Although simpler methods would have given the same result in this case, the 
method is very powerful, since it may be used on a solution with no further reduction 
from the contour integral form. 

6. Semi-infinite slab, 7>,, polynomial in time, initial temperature zero. Let us 
consider the following equation 


C(2n)! a+ 2cosz\*. 
T2.¢ = (— 1)* f ( ) eterno, (6.1) 
201 P a+2 


An analysis similar to that in §5 shows that 7,,9=0 and also that the integrand 


vanishes over the paths M and UN. 

When x=0, the integrand is even-valued in z, and therefore has the same value 
over path Q as over P. To,; may then be equated to half the sum of the residues at the 
poles enclosed by paths P and Q. From (12.8) this becomes 


To,e = Ean(t), (6.2) 


where £2,(¢) is a polynomial in ¢ of the mth degree. From proper combinations of these 
polynomials, contour integrals are readily derived for problems in which the boundary 


temperatures are arbitrary polynomials in time. 
For the particular case where the boundary temperature is linear in time and 
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therefore given by To,,= C&(t) = Ct/(a+2), as follows from (12.10), the contour solu- 
tion is given by, 





. a+2cosz\‘_ dz 
T:24= —-— ——— ] e** — (6.3) 
a) a+2 z* 
We indent the origin, integrate as in §5 to obtain 
2C a a 2cosw\! . dw C */a + 2 cos ee*\! exp [ixee’*] 
a ——— } sin wz — —- — : do 
wd. a she y w? ro a+2 €*e7*¢ 


The first integral, after the expansion of (a+2 cos w)‘ in terms of the polynomials and 
regrouping as in (5.2), becomes 


Fr = 
sin w(x + 7) P(t) — - 
” we oa ane nl e w* 
By integrating by parts and keeping all terms which do not drop out as e—0, the con- 
tribution of this integral is found to be 


1 mw 208 C x tor wf etek , 
SE a st ree(l sm Ws 7) —* 
Te m(a+ 2)? woe iia is ‘ w 


The second integral is expanded in terms of ¢€ and ¢, and all terms are retained which 
do not drop out at the integration limits or as e—0. We then have for this integral 


Cz CH Ct 
Te 2 a+2 








,>- 
The sum of J; and J, with e—0 then gives for the solution 


t x? 1 r=+t 
Te4=C + — —- ———_ atr)Pi(d)ixtr | (6.4) 

l=; 2 2%a+2)! 2¢ +e) 
where the symbols {x+r} have the same meaning as in (5.3). 

An analysis similar to that of §5 shows that if a20, the numerical solution ap- 
proaches Ct/(a+2) asymptotically as t— ©. This result also follows from the analytic 
theory. 

7. Finite slab, length /, boundaries x = 0 and x=/ held at constant temperatures 
uy and u; respectively, initial temperature zero. Let us consider the equation 

1 a + 2 cos z\' . dz 
T24=— (== =*) [A(z) cos zx + B(z) sin zx]—> (7.1) 

rt a+2 z 
where the functions A(z) and B(z) must be determined to make (7.1) satisfy the 
boundary conditions T»,,=%o and 7),,=«;. Referring to (5.1) and placing x =0, we 


see that 

u a 2 cos z\' dz 

s — =) __ (7.2) 
rt a+2 Zz 


Therefore if A(z) =u and B(z) =(u;—wo cos 2/)/sin zl, (7.1) reduces to u» at x =0 and 
to u; at x =/. The solution is therefore 
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1 a+ 2 cos 2\' ui— Ugcoszl . dz 
T 24 = —f —_——_—_———_ %6o COS 3x -- ———-. sin 8x ]— (7.3) 
P a+2 sin 2l Z 


The integrand of (7.3) is an even-valued function of z and therefore the solution 
may be equated to one-half the sum of the residues at poles enclosed by paths P 
and Q (there are no contributions from the paths M and JN). Poles occur at z=n7/I, 
n any integer. After evaluating the residues and simplifying the result, we have 


x 
T 2.4 = Uo t+ (ui — Uo) 7 





ja ( LD) fe 2 /l 
+ — ¥> (— 1)"(ur — uo cos nm) = cota “(== oo ey ; (7.4) 
WT 1,2 n a+2 


To express this solution in terms of the polynomials, it is necessary to expand 


(a+2 cos (nm) /l)* as in (12.7), 


x 
T 2,2 = Uo + (uz — Uo) ri 


4 > (— 1)% sin (nmx/ Dd [?. () + 2P,a(2) nT 
— — 1)"(u; — 1 oe cos — 
ta + 2)! as ui Uo COS NT - t—1 Ss i 





nit 
+ +++ + 2Po(t) cos =| 


Combining the trigonometric terms and making use of the identity P ._,(t)=P.4,(¢), 
we obtain the solution in the form 


x 


T 2,4 = Uo + (ui — Uo) " 


r=+t 
a > Pus) ed (u%, — Uo COS mm) sin un = . 
am(a+2)* -; n 
However from the initial condition 7,,9=0, it follows that at ¢=0, the resulting 
sine expansion of (7.4) must be equal to —u 9—(u;—uo)(x/l). Therefore the infinite 
series in the double summation of (7.5) must equal — F(x+r), where F(x) is the peri- 
odic sine expansion of u9+(u;—4u9)(x/l). The solution then becomes 








(7.5) 


1 r=+t 


x 
T 2.4 = tot (uy > uo) « ial (a+ 2) zu P, (dF (x + r). (7.6) 


Figure 3 shows the graph of the periodic sine expansion of F(x) =uo+(uz—o) (x/I) 
which applies to (7.6). In applying (7.6) it is 
usually simpler to plot the function and then 
pick off the different values of F(x) required 
in the summation. It will be noticed that 
when x+r is a multiple of J, the value of 

F(x+r) is zero. 
For finite slab problems, it is usually 
Fie. 3. possible to get a solution in the form of a 
finite Fourier series, in addition to the poly- 





’ 
[a 
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nomial expansion. The derivation of such series is not difficult, but is too long to be 
given here in full. The finite series for §7 is 
3 na x (° + 2 cos my 


T + ( yt > A, si (7.7) 
zt = —%)—— . sa — ; : 
” a ™ l 1 . l a+2 





where 
sin (nw/l)(uzy — uo COS mm) 


l(1 — cos (nx/1)) 


The general method used for obtaining finite series, such as (7.7), consists of a 
replacement of the infinite Fourier series as in (7.4) by a finite number of. terms of 
the same type, such that at t=0 the finite series reduces to the same function as given 
by the infinite series at ¢=0. Such expansions are possible since in numerical methods 
the initial temperature must be specified only at a finite number of points. The coeffi- 
cients for the terms of such finite series are given by® 


(7.8) 





n 


ro, nex 
A =—)> F(x) sin 


z= 





’ 


where F(x) is the function at t=0 over the interval 0 to /. The A, of (7.8) were cal- 
culated as above with F(x) =uo+(u,;—o)(x/l). 

From an inspection of (7.7) we see that as t+ the solution converges provided 
| {a+2 cos (mr/l)}/(a +2)| <1. A simple analysis yields the following criterion for 
convergence 


1) a+2 2 2cos?(x/2l) 1 even, 


(7.9) 
2) a+222 cos? (x/l) 1 odd. 


Equation (7.7) also shows that the numerical solution approaches uo + (u; — uo) (x/l) 
as t+, which is the analytic steady-state solution for the same boundary conditions. 

8. Finite slab, length /, insulated at x=0, held at constant temperature 1; at 
x =/, initial temperature zero. The boundary conditions are 


Ti = ui, (8.1) 


and from Eq. (2.7), 
2T 1 4- aT 9 + 
— 10-1 + OT t-1 (8.2) 
a+2 





Imposing these conditions on the general contour integral (4.2), we find that 
A(z) =u;/cos zl and B(z) =0. The solution is therefore 


U1 a+ 2cosz2\‘ cossx dz 
T 23 = — ( ) ae (8.3) 
a+2 


Tid p coszl 2 
Poles occur at z=0 and z=(2n+1)2/2l. Evaluating the integral as in §7 and regroup- 
ing terms, we obtain the solution 





5W. E. Byerly, Fourier’s series and spherical, cylindrical and ellipsoidal harmonics, Ginn Co., 
Boston, 1895, pp. 30-35. 
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ye uf n 4 > yy yontnn(* + 2. cos (nx/2l) ) cos are, (8.4) 


T 1,3 ao+2 n 
An analysis similar to that of §7 yields the polynomial expansion 
1 x 
T 34 = ak _- ———_ Pri lDF (x + | (8.5) 
(a + 2)° r=—t 


where F(x) is the cosine Fourier expansion of unity from 0 to/ and minus unity from 
1 to 21, as is evident from (8.4) and the initial condition 7,,9=0. It will be noted that 
the value of F(x+r) is zero when x++7 is an odd multiple of /. 

The finite Fourier series solution is found to be 


35-1 nx(x — I) (- + 2 cos (nx/2l) | 
T ;=> 1 An i ’ 8.6 
nf + X sin rr ee ) (8.6) 





where 
sin (n2/2l) 
A, =0 neven, A, = ——_._._ nn Od. (8.7) 
l(1 — cos (nw/2l)) 


From (8.6) it follows that the convergence criterion is: 


1) a+2 2 2 cos? (4/4) 1 even, 


2) a+2 2 2 cos? (x/2l) 1 odd. (8.8) 


Also from (8.6) it follows that as t— © , the numerical steady-state solution becomes u, 
which is also the analytic steady-state solution for the same problem. 

9. Finite slab, length /, constant energy input qg at x=/, temperature kept at 
zero at x =0, initial temperature zero. The boundary conditions are (from Eq. (2.5)) 


2T y~-1,4-1 + aT y 1-1 + 22 
>] 
a+2 
where Q=gAx/k. Imposing these conditions on (4.2) we find that A(z)=0 and 
B(z) =Q/(sin z cos 2/). The solution is therefore 
sin 2x a + 2 cos z\' dz 
@ iceteuee. a 
a+2 z 
Poles occur at z=0 and z=(2n+1)z/2/1 (the set z=nz, n 0 does not constitute 
poles, as the terms sin zx in the numerator also vanish at these points since x is re- 
stricted to integral values). After evaluating the residues, regrouping terms and sim- 
plifying, we see that the solution becomes 
sin (nrx/2l) (- + 2 cos i 
a+2 , 


(9.1) To, = 0, (9.2) 





Ti = 





T24 = — . 
mid p sin z cos 2l 





T 24 = o| « + £ > (-— 1) (er) /3 (9.4) 


T 13 n sin (nx/21) 


The polynomial expansion is found to be, 


1 r=+t 
= —_— ————_ t 9.5 
Tea =O[1-— E Palorte +], (9.5) 
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where F(x) is the sine expansion of x from 0 to / and of 2/—x from / to 21, as follows 
from (9.4) and the initial condition 7,,5=0. 
The finite Fourier series is 











a nex {a+ 2 cos (nr/2l)\! 
T24 = o| = — >> A, sin ( ov ‘y'| (9.6) 
1 21 a + 2 
where 
(— 1)@-pr2 
A, =0O wn even, A, = n odd. (9.7) 
l(1 — cos (nw/2l)) 
From (9.6) the convergence criterion is found to be 
Tr 
1) a+2 2 2 cos?— / even, 
4] 
Tv 
2) s+ 2eiot— 1 odd. (9.8) 


From (9.6) it also follows that as > the numerical solution approaches Qx, which 
is equivalent to the analytic steady-state solution for the same boundary conditions. 

10. Convection at a boundary. When the condition (2.3) is imposed on the gen- 
eral contour integral (4.2), A(z) and B(z) are generally such that the evaluation of 
the resulting integrals is difficult, due either to the uncertain nature of the poles, or 
to the evaluation of a complicated infinite integral. 

For the semi-infinite slab with convection into temperature T,, transfer coeffi- 
cient h, the boundary condition is given by (2.3) 





98.0.0 4 le: — Din td 107 
a il eT ns. ° iv hath (10.1) 
a+2 
We assume a contour integral solution of the form, 
1 a+2cosz\!_ dz 
T2.2=— | A(z) (6) e#7 —. (10.2) 
Tid p a+2 z 


Imposing the condition (10.1) on (10.2) we see that A(z) = NT.,/(N—1 sin 2), and the 
solution may therefore be written, 


Ni« a+ 2 cos z\¢ ett dz 
T34= f ) —- (10.3) 
P 


rt a+2 N—isinz z 


The integral (10.3) is to be evaluated in terms of an integral along the real axis 
indented at the origin, plus the sum of the residues at the poles enclosed by the path 
P and the real axis. Aside from the root z=0, the denominator includes roots from 
the term N—i sin z, which are found to be z= —i log (\/N?4+1+N) and the infinite 
set = —i log (\/N?4+1—N)+(2n+1)z, n any integer. However, in the loop of inte- 
eration considered, only the residues at the poles = —i log (1/N?+1—N)+(2n+1)r 
are evaluated, since the other poles lie outside the loop. (N is taken greater than zero, 
otherwise the boundary would be insulated and no heat would flow, giving the trivial 


solution 7,,,=0.) 
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Following the analysis of (5.1), integrating along the indented real axis, and then 


letting e—0, we obtain the solution in the form 


r r E 2N C + 2 cos “\(- sin wx + sin w cos =) dw 
: 0 a+2 N? + sin? w w 





T 
+ 2ni >, residues | (10.4) 


By a somewhat tedious but straightforward analysis, the residue term may be evalu- 
ated and simplified to yield, 





4NT, lo ‘N?+1-—-—N) /a—2VN?+1 
dni D; Res. = 6 (N ( 2 


(/N?+ 1 — N)? 
VN? +1 a+ 2 =) 


COS NTX 


i ee -+ (10.5) 


1,3 nr? + log? (\/N? + 1 — N) 





By choosing a combination of known Fourier expansions for the hyperbolic sine and 
cosine, with further reduction, and recalling that in numerical analysis x is always an 
integer, we can write (10.5) in the form 


—2/N?4+1 /N?4+1—-1 
ani D Res. = - T.(* a ) ew = VIFF +(- nd ). 
a+2 VN? + 1 











The final solution therefore becomes 


r , E = f ic + 2cos =\(- N sin wx + sin w cos =*) dw 
pois rT Jo a re 2 N? + sin? w w 


wei — 2/N?+1\! gerbe Fike 
-(~ ; ——) es yaw — S/N? + =|. (10.6) 
/N? + 1 a+2 


Without evaluating the infinite integral (which is convergent when t-~ provided 
a=0) we see by inspection of the factor [(a—2./N?+1)/(a+2)]* that the criterion 
for convergence as t— & is, 











a = J/N?4+ 1-1. (10.7) 


It will be noted that Schmidt’s equation, where a=0, will not yield convergent 
answers if the boundary expression (2.3) is used, since by (10.7) @ cannot be zero. 
It also follows from (10.6) that the numerical steady-state solution becomes 7,, which 
is also the solution predicted analytically for the same problem. 

A polynomial expansion may be obtained from (10.6) by using (12.7) to expand 
the term (a+2 cos w)‘. However, the analysis required to evaluate the resulting defi- 
nite integrals is so involved that it is not worthwhile to include the expansion here. 

Contour integral solutions are readily set up for finite slab problems involving 
convection at either or both boundaries, but it is generally difficult to evaluate these 
integrals. However, without evaluating the integrals, but using the method outlined 
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in (5.3), we can show that if the numerical work converges at all, it converges to the 
correct steady-state distribution as t>. 

11. Convergence to analytic solution. The theory developed so far has been con- 
cerned only with numerical solutions as actually applied. The convergence problems 
already discussed have shown what values of a are necessary to insure non-divergent 
numerical answers, and also that as the time increases, numerical solutions approach 
the same steady-state distribution as that given analytically for the same problem. 
These results hold when the space and time increments Ax and A? are finite, non-zero 
quantities. 

It is also possible to show that as the arbitrary increments Ax and At become very 
small, the numerical solutions approach the true analytic solution for all values of x 
and #, and attain this solution in the limit. The formal procedure used to demonstrate 
this limiting convergence consists of a demonstration that the contour integrals de- 
rived for the numerical solutions transform into already known contour integral solu- 
tions for the corresponding analytic treatment, as Ax and At approach zero. Formal 
proofs of this convergence will be given for three of the numerical examples already 
discussed. The proofs for other examples are very similar to the ones given here. 

The three examples considered with their analytic contour integral solutions are :* 

1) Semi-infinite slab, end x’=0 kept at temperature uo, initial temperature zero, 





uo a , dz 
_ T(x’, tf) =— |] ette’e- et? —. (11.1) 
Tid p 2 
2) Semi-infinite slab, convection at x’=0 into a medium of constant temperature 
T, initial temperature zero, 
hT. ett’ east’ dz’ 
T(x’, ¢’) = - f n—areemmeeermeararty mtomee © (11.2) 
kni Jp (h/k) — iz 


3) Finite slab, length 1’, end x’=0 held at zero temperature, constant energy 
input g at x’=/’, initial temperature zero, 
q sin 2x’ dz 
Td, if) ee beer’ (11.3) 
kriJd p cos zl’ 3? 

Equation (11.3) is not given in Ref. 6, but we can easily derive it by imposing the 
boundary condition 07/dx —q/k =0 on the general contour integral considered there. 

The notation used in the above equations differs from that used in the present 
paper. In order to express these equations in our notation, it is necessary that the 
complex variable a be replaced by z, thermal diffusivity K by a, and convection co- 
efficient h, by h/k. 

T(x’, t’) has been used to denote the analytic solution at the point x’ and the 
time t’, as distinguished from 7,,:, the numerical solution at point x (in units of Ax) 
and time ¢ (in units of At). The path P is the same for both numerical and analytic 
solutions, and is the limiting path allowable for the analytic integrals. 

In the application of numerical solutions, position and time variables as well as 
slab lengths are expressed in terms of the arbitrary time and space increments At 
and Ax. In order to discuss the convergence to analytic solutions, it is necessary to 
express these quantities in terms of the absolute units used analytically. If M arbi- 


6 H. S. Carsiaw, The conduction of heat, Macmillan Co., New York, ed. 2, 1921, pp. 97-99. 
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trary space increments equal one absolute unit, then the absolute space position x’ 
is given by, 
zs’ =k/M or «= Mz’. (11.4) 
From the relationship At =(Ax)?/(a+2)a@ of (2.2), it follows that the time value in 
absolute units is 
= 1/M%{a+2)a or t= Mat 2at’. (11.5) 

To show that the numerical solution approaches the analytic solution for Case (1), 
when Ax and At approach zero, we substitute (11.4) and (11.5) into the numerical 
solution (5.1), obtaining 





u a + 2 cos 2\¥?2(et2) at! dz 
= ( ; = Praia (11.6) 
Tid p a+2 z 
By replacing the variable z by 2’/M, we can write (11.6) in the form 
u a+ 2 cos (2//M)\MPet2at dz’ s. 
ae ( ) gis’? (11.7) 
Tid p a+2 


The factor [(a+2 cos 2’/M)/(a+2) ]¥*(s+%« may be approximated by [1—(z’)? 
/(a+2)M?]¥*(2+2)<" when M is large, and in the limit becomes exp[—(z’)*at’ ]. The 
limit of the numerical expression (11.7) as Ax and At approach zero (or as M— ~~) 


therefore becomes, 
/ 


uo ee , , dz 
er eee Ys ett'2' gale Pe any (11 . 8) 


lim Az 0 Tid p 2! 


which is identical with the analytic solution (11.1). 
With convection at a boundary or constant energy input, the values N =hAx/k 


of (2.4) and Q =qAx/k of (2.6) become, 


h 
aes (11.9) oer (11.10) 


kM kM 
In Case (2) after substituting (11.4), (11.5) and (11.9) into the numerical solution 
(10.3), and then replacing the variable z by z’/M, we obtain 


hT, KG + 2 cos a ef’s’ dz’ (11.11) 
 kwiM Jp a+2 (h/kM), — i sin (2//M) 2’ 
As before, the term [(a+2 cos 2’/M)/(a+2) ]¥*(¢+2«" becomes exp[— (z’)?at’ | in the 
limit. The term 7 sin(z’/M) may be approximated by iz’/M when M is large, We 
make these changes, cancel the M outside the integral with those in the term 
(h/kM) —(L2’/M), and let M approach infinity; then (11.11) becomes 


ET. ett’ 2’ eG’ tat! dz’ 
a — f Teeedcaees eta? (11.12) 
lim Az—0 kri Jp (h/k) — iz’ 2’ 











z,t 


which is the analytic expression (11.2). 
In Case (3), after substituting in (9.3) from (11.4), (11.5), (11.10), introducing 
the absolute length /’ given by /= MI’, and changing the variable z to z2’/ M, we obtain 


a 2¢ 1/Yf)\M2(at2)at? de! 
g f sin 2'x (- + 2 cos (2’/ ‘) & (11.13) 
F 


et kwiM J p sin (z’/M) cos 2'I! a+2 








Z 
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In the limit, the term sin(z’/M) may be replaced by 2’/M and as in (11.11), the solu- 


tion may be written in the form 
q sin 2’ x’ F 
lw eo ———  ¢- (#’ Rat! , 
lim Az—0 krid p cos 2'l’ (z’)? 


(11.14) 





which is identical with the analytic solution (11.3). 

It will be noted that the restrictions on a as given by (9.8) become a20 when / 
(the number of units of Ax in the slab) becomes large. Hence the proof for this con-: 
vergence to the analytic solution holds only when a20. An analysis of the other ex- 
amples treated in this paper shows that for this limiting convergence, all criteria re- 
duce toa20. 

12. Appendix. Properties of the polynomials P,(t). The polynomials P,(t) are 
defined as the coefficients of z* in the expansion of the trinomial (1+az+27)*, and are 
therefore functions of r, t and the modulus a. Two identities follow readily, the first 
by definition, the second by setting z=1 in (12.1), 


r=+t 


D> Prcé)stt? = (1 + az + 2%)', (12.1) 
r=+t 
> Put) = (a+ 2)4. (12.2) 
r=—t 


By expanding the trinomial in the form [(1+az)+2*]‘ and collecting coefficients, 
we obtain an explicit formula for P,(¢), 


rav-()(Jer (rere (DG) eas 


The polynomials may be expressed as definite integrals in the following way. By 
definition and Cauchy’s theorem 


1 », 22 
Pt) = <= f Ot ost 


where C is a simple closed contour about the origin. By the choice of C as a circle of 
unit radius, center at the origin, it follows that z=e'* and 


1 Qn 
Pi) = —{ (a + 2 cos ¢)' cos (t — r)ddd. (12.4) 

rd 9 
Writing the equality (1-++az+2s*)'=(1+az+2?)'1(1+az+2?) and collecting coeff- 
cients of z* on both sides of the equation, we have the following recursion formula, 
which may be used for rapid calculation of the polynomials, 


P(t) = PAt — 1) + aP,(t — 1) + P,-2(t — 1). (12.5) 

From (12.3), (12.4) and (12.5), it follows that 
P,(é) = Pi(d), (12.6a) 
P,(0) = 1, P,(é) = 1, P,(0) = 0, r #0. (12. 6b) 


As an example, the polynomials up to t=3 have been worked out for modulus 
a=3 using (12.5) and (12.6), and are shown in Table 1. Thus, P,(3)=30 and P2(2) 


=11. 
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TABLE 1. 
t r | 0 1 2 3 | 4 | 5 6 
0 ! | | 
1 1 3 1 | 
2 i 6 11 | 6 | 1 
3 1 9 30 45 30 9 1 


To construct a polynomial array, we start with P»(0) =1. The polynomials follow- 
ing are calculated successively by use of the recursion formula (12.5). As a specific 
example, from the formula P3(3) = P3(2)+aP2(2)+P:1(2), we have on substituting the 
values presumably already calculated for t=2, P3(3) =6+3X11+6=45. 

An important identity may be established as follows: we let (a+2 cos 6)‘ 
=)-!_ A, cos 70. From (12.4), 


2r t 
P,(t) = f - A,, cos n@ cos (t — r)ddd, 
0 


n=0 


from which it follows that A ;_,=2P,(t), r¥#t, and Ap =P,(t). Therefore 
(a + 2 cos 6) = P,(t) + 2P,_;(t) cos0 +--+ + 2P:_,(t) cosr6+--:- 








+ 2P (t) cos 2. (42.7) 
The polynomials £,(¢) are defined by, 
(— 1)" (2n)! a+2cosz\‘ dz 
fon(#) = fi ) , (12.8) 
2 2mi Joc a+2 giott 


where C is a simple closed contour about the origin. By Cauchy’s theorem, and evalu- 
ation of the residue at the pole z=0, (12.8) becomes, 


(— 1)*f d* fa + 2 cos @ | 
n(é) = ——_——————— A 12.9 
tnt) =§ [a ) | (12.9) 


The first three polynomials, evaluated from (12.9) are, 








t ) 
(0) = ( (a +2), | 


E(t) = (—) (a+ 2)-1+ 12( Jc + 2)-2, 


(t) = ( ) (+2714 00 ( a +224 360( ) (a + 2)-4 


| 
+ = (12.10) 
| 

| 
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—NOTES— 


ON PLANE ELASTIC STRAIN IN DOUBLY-CONNECTED DOMAINS* 
By W. PRAGER (Brown University) 


1. Introduction. The stresses associated with a state of plane elastic strain can 
be expressed in terms of the second derivatives of Airy’s stress function. If x1, x2, x3 
are rectangular Cartesian coordinates, the axis of x3 being normal to the plane 
of strain, the stress function $(x:, x2) satisfies the differential equation A*¢#=0 
(A =0?/dx?+0?/dy?), and the given stresses on the boundary determine the tangent 
planes of the stress surface x3=@(x1, x2) at all points of the boundary, when one such 
tangent plane is known for each bounding curve. In the case of a singly-connected 
domain only one such tangent plane must be known, and it can be chosen arbitrarily 
because the stresses define the stress function only to within an arbitrary linear func- 
tion of x; and x2. In the case of a doubly-connected domain, however, two such tan- 
gent planes must be known, and only one of them can be chosen arbitrarily. This 
paper is concerned with the determination of the second tangent plane in the case 
where one of the boundary curves is free from loads. Equations from which this 
tangent plane can be determined, were derived by J. H. Michell from the condition 
that the displacements must be single-valued. In the present paper it will be shown 
that Michell’s equations are the natural boundary conditions of the variational prob- 
lem for the stress function. This remark is of importance when the direct methods of 
the calculus of variations are used to determine the stress function for a doubly- 
connected domain.? 

2. Notations. Basic relations. Throughout this paper Latin subscripts will have 
the range 1, 2, 3, Greek subscripts the range 1, 2, and the summation convention 
for repeated subscripts will be used. The rectangular Cartesian coordinates x; are 
chosen so that the axis of x; is normal to the plane of strain; the position of the origin 
and the directions of x; and x2 are arbitrary. Let e;; be the strain tensor and s;; the 
reduced stress tensor, i.e. the stress tensor divided by Young’s modulus. The stress- 
strain relations can then be written in the form 


6g = (1 + o) 55 — oSKidi 5, (1) 


where o denotes Poisson’s ratio, and 6;; is the Kronecker delta. For the state of plane 
strain under consideration the condition that e33;=0 gives 


533 = OS yy. (2) 
The equations of equilibrium in the plane of strain are 


* Received Aug. 7, 1945. 

1 J. H. Michell, Proc. London Math. Soc. (1) 31, 100-146 (1899), Eqs. (13). 

2 The necessity of investigating the relations between the natural boundary conditions and Michell’s 
equations arose in connection with work done under a contract in Applied Mechanics for Watertown 
Arsenal. The author is indebted to the authorities of Watertown Arsenal for the release of this note for 
publication. 
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Sya,y = 0, (3) 


where the comma followed by the subscript y denotes partial differentiation with re- 
spect to x,. Equation (3) can be satisfied by setting 


Sap = Ear€ uP rus (4) 


where ¢=9(x1, x2) is Airy’s stress function, and €1=€2=0, é2=—é€:=1. Since 
€ad€ay, = 5y,, the invariant Seq equals A®=@.ea. For a state of plane elastic strain 
ASaa =0, or A*d =0, i.e. the stress function is biharmonic. 

On the boundary of the domain under consideration the surface stresses f, are 
given. If m, is the unit vector along the outward normal of the boundary, we have 
Sa =SyaNy = €yr€aph,ruNy. Now €,,n,=t,, the unit vector of the tangent of the bound- 
ary. Accordingly, 
Soe = €apP uh = €an0P,,/9S, (5) 
where 0/ds denotes differentiation in the direction of the tangent vector ¢,. Multiply- 
ing both sides of (5) by €ag and integrating along the boundary, we obtain 


6918) = eas f fals)ds + 4.010). (6) 


The given surface forces f, are thus seen to determine the gradient ¢,,(s) of the stress 
function along a bounding curve, when the gradient ¢,(0) at one point of this curve 
is known. In other terms, the stress function @ and its normal derivative 0¢/0n are 
defined at all points of a bounding curve, when ¢ and its gradient are known at a 
single point of this curve. 

If, in particular, one of the bounding curves is free from loads, the stress function ¢ 
and its normal derivative along this curve equal the values of a linear function dax2+b 
and of its normal derivative. Establishing the boundary conditions for the stress 
function along a bounding curve which is free from loads is therefore equivalent to 
determining the three coefficients a, a2, b of this linear function. 

3. The variational problem for the stress function. To the strain e;; and the re- 
duced stress s;; corresponds the reduced elastic energy U = }e;;5;;. In the case of plane 
strain this energy equals 


U = $3[(1 + o)sijsi3 — oss 5;] = 4(1 + 0) [SapSap — oSaaSpe, (7) 


in view of Eqs. (1) and (2). In terms of the stress function introduced in (4) the en- 
ergy is 

U= 3(1 + a) Cee onal od aah pp]. (8) 
According to the variational principle for the stresses the stress function correspond- 
ing to certain boundary conditions is then singled out from amongst all functions 
which fulfill these boundary conditions and admit continuous derivatives up to the 
fourth order, by the fact that it minimizes the integral 


V= f [¢.096,08 — 0,00 ,68 \dw, (9) 


where dw denotes the element of area, and the integration is extended over the entire 
domain. The condition 6V =0 leads to 
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f [¢ 295, ap = 59 aad ap |deo = 0 (10) 


or 


ts 0) | baassdédes -(i- 0) fo assbneis 


+ f $.0086.0meds —- «| 1986 aed = (). (11) 


In the case of a doubly-connected domain with loads on one bounding curve only, 
the stress function ¢ and its gradient can be considered as given on the loaded bound- 
ing curve. This curve does therefore not furnish any contribution to the line integrals 
in (11). On the other bounding curve, we have ¢=a@,%2+6 and $,4=d,. In addition 
to the differential equation for the stress function, ¢,2es3 =0 or A*@ =0, Eq. (11) thus 
gives the following equation which must be fulfilled on the load-free boundary: 


say] ( _ 0) f .asttads = frames + of &o0nds 
+ x| (1 - 0) fo atonads = 0. (12) 


Since 5a, and 6b are independent, the expressions in brackets must vanish separately. 
The second integral in the first bracket can be transformed as follows 


J srmads = f bseyde = f osama 


The first bracket can therefore be written as 


(1 ae o)| ff a0r2neds — f eaonras | 


With the use of ”,=€,ata, the second integral can be further transformed as follows 


J ¢.s0mds cou f d.sntds = conf band 


_ cou f d.ontads = ew f @ ppXateds. 


Il 


Equation (12) is thus equivalent to 
rs) re] 
J | #1 = Go) + erate = (49) Jas = 0 (13) 
on Os 


and 


fz (Ad)ds = 0. (14) 
on 


The scalar equivalents of (13) are 
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a 


= — 
f as — (46) + — We) bis'= 0 (15) 
On Os e 


and 








r 9 3 ae 
f | #2 — (A) — a= (a9) [as = 0. (16) 
Equations (14), (15) and (16) are Michell’s conditions which are thus seen to be the 
natural boundary conditions of the variational problem for the stress function. The 
manner in which these equations are used in determining ¢ is obvious. Let ¢o, $1, $2, 3 
be the biharmonic functions defined by the following boundary conditions: 
1) do and 0¢0/dn have the prescribed boundary values on the loaded boundary 
curve C, and vanish on the other boundary curve C2; 
2) $1 =0¢g1/On =0 on Ci, 
¢$1=%, and 0¢;/0n =n, on (2; 
3) ob: =0¢2/dn =0 on Gi, 
doo = x2 and 0¢2/dn = nz on C2; 
4) d3=0¢;/0n=0 on Ci, 
3=1 and 0¢,/dn =0 on Cy. 
Substituting 
d = go + A191 + A2h2 + Dds 


into Eqs. (14), (15) and (16), we obtain three linear equations from which a, a2 and 


b can be determined. 


THE CAPACITY OF TWIN CABLE—II* 
By J. W. CRAGGS anp C, J. TRANTER (Military College of Science, Stoke-on-Trent, England) 


1, Introduction. In a recent paper! (subsequently referred to as “I”) we have 
given a method for determining the capacity of two circular wires surrounded by con- 
centric touching dielectric sheaths. The present note gives the extension of the method 
to the case in which the dielectric sheaths are not in contact. The problem considered 
is the symmetrical one of two infinite parallel circular wires each of radius R; sur- 
rounded by concentric sheaths of radius Re and dielectric constant K,, the distance 
between the centers of the wires being 2L(L>R-,). The dielectric constant of the sur- 
rounding medium is taken as Ke. 

2. The equations for solution. In line with the treatment in “I” we replace R. 
by unity, Ri/R2 by a and L/R; by s; we also write Ki/K2=K. The potentials Vi, V2 
must therefore satisfy (i) the differential equations 

V*V; = 0, esrs t, (1) 
V*V2, = 0, r = 1, x2 0, (2) 


and (2) the boundary conditions 
V; = 1, (3) 


* Received June 19, 1945. 
1 J. W. Craggsand C. J. Tranter, The capacity of twin cable, Quart. Appl. Math. 3, 268-272 (1945). 
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when 7=a, 
V: = Vs,, (4) KaV;/dr = 8V2/dr, (5) 


when r=1, 
Ve = 0, (6) 


when x =0. Here V? is Laplace’s operator in two dimensions and the coordinate sys- 


tems are as shown in Fig. 1. 


(_ 
Nut 











Fie. 1. 


3. The analytical solution. As in “I” we write 


r 2) rT n a n 
Vi; =1+ Blog— + y {(<) ~ (=) bo, cos 6. (7) 
a n=] a T 


The conformal transformation fer the region r>1, x >0 can be written 


re? + eH 








£ — in = log — (8) 
. re + e-# 
where 
uw = log (s + Vs? — 1). (9) 
The boundaries r=1, x =0 then become =y, £=0 respectively. 
Since V2 is odd in & and even and periodic in 7, we write 
V2 = DE + > dm sinh mé cos mn. (10) 


m=1 


The constants B, b, of (7) and D, d,, of (10) are now to be determined from the 


boundary conditions (4) and (5). 
On the boundary r=1(£=y), we find from (8) and (9) 


1 + cosh pu cos 6 (11) 


cos y = ——-———_ 
cos 6 + cosh u 
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so that 0 <@ <7 corresponds to 0 =n S7, and 

















OV 0& OV On OV — sinh pu OV 12) 
eect wer cca: Seles ae paaeee (12 
or Or o€& 06 dé cos@-+ coshy dé 
Thus (4) and (5) give 
a) 1 — q?" a) , 
1— Bloga+ >> ———b, cos n0 = Du + >> dm sinh mp cos mn, (13) 
n=1 a* m=] 
© 1+ a 
KB+K)>, — nb, cos n@ 
n=1 a” 
— sinh yu oe 
= {p + >° mdm cosh mu cos mn) ‘ (14) 
cos 6 + cosh yu m=1 
Multiplying (13) by cos mn(m=0, 1, 2, - - - ) and integrating with respect to 7 
from 0 to 7, we have 
Oo 1 = q?* ; 
Du =1-— Bloga+t+ z, (— 1)"e-™ —— },, (15) 
a a” 
since 
f cos nOdn = (— 1)"re-™, 
0 
and 
 -) 1 _— 2n 
dm sinh mp = > e~"# — bp,Im(n), (16) 
n=1 a” 
where 
2 . ” 
I,(n) = — ew f cos n@ cos mndn. (17) 
T 0 
Similar treatment of (14) gives for B, bn 
KB=-D (18) 
and 
_1+a™* = 
K — -nb, = 2(— 1)"*e-™D — em >> mdm cosh mpl m(n). (19) 
a” m=1 
Expansion of cos m7 in (17) in terms of u=(1+e—**+2e-* cos 6)—' leads to 
Im(N) = (ae 1) mtn) (— 1)? a n+P-1C ¢(n—2p)u, (20) 
p=0 
Eliminating D, b, from equations (15), (16), (18) and (19) we have 
1 oO 
B log a — 1 — KBu = 2BS + — > mdnam cosh mp, (21) 
_ m=) 
. 1 ~ 
— d,sinh pp = Ba, + — >> md mA mp COSh mp, (22) 


m=1 
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where 
2/1 — a*"\ e 2 
S = x ( ); , | 
n=1 1 + re n 
«© 1 — a2" e727 
= 2 — 1)*+1| ————__} 7 , 23 
ay = 23, (— (rm, (23) 
| 
i) 1 =s 2n e727" 
Am, = ——— J T,,(n)I . 
=X (2) mmc) — | 


Following the procedure of “J” we retain only a finite number p of the coefficients 
d». Writing 


, 


Ym = — tanh mu, (24) 
m 


and eliminating md,, cosh mp between equations (21), (22) we find 








| Auty Ai it Ai, 1 

| Ag Ast ¥2°** Aop a2 

= 0. (25) 
A p1 A p2 s2° A pp t+ Yp ap 

1 

ay a2 tee ap 2(— — tog «+ Ku) + 45 


The capacity is then given by —}K,B. 

4. Alternative method of solution. The above treatment provides a satisfactory 
basis of computation when K 21. For completeness it is interesting to notice that, 
when K <1, more rapid convergence to the true solution is obtained by eliminating 
D and d,, from equations (13), (14) by treating (13) as a Fourier series in @ and (14) 
as one in 7. 


ON A. A. POPOFF’S METHOD OF INTEGRATION BY 
MEANS OF ORTHOGONALITY FOCI* 


By HOWARD A. ROBINSON (Research Laboratories, Armstrong Cork Company) 


In a recently published paper! a method is given which allows a marked reduction 
of the work necessary in computing the tristimulus values necessary in color specifi- 
cation work. The three tristimulus values are defined by the following relations: 


X = [ £:0)20)RO)m, Y = f 2.0)s=)RO)A, Z = f E.0)2-)RO)A, 


where £,(A) are tabulated relative energy functions of a known light source L, #(A), 
j(A), 2(A) are tabulated luminosity functions and R(A) are the experimentally meas- 


* Received August 9, 1945 
1 Quart. Appl. Math., 3, 166-174 (1945). 








384 BOOK REVIEWS [Vol. III, No. 4 


ured curves of reflecting power or percent transmission as a function of the wave 
length. The product functions E,%, E,j and E,2 have also been tabulated and may 
be considered as the ¢;(x) in the original article. Work is now under way in setting up 
the necessary scales for the colormetric computations and will be published elsewhere. 

It is pointed out in the introduction to Popoff’s paper that the method requires 
the construction of certain diagrams, called scales, showing the abscissas of the cen- 
troids of certain areas associated with ¢;(x). Thus, operation (b) in Section 2 contains 
some unnecessary work, since it is unnecessary to find the centroids 4,, 5,, - - - , only 
the abscissas of these centroids being required. 


BOOK REVIEWS 


Elementary eleciric-circuit theory. By Richard H. Frazier. McGraw-Hill Book Com- 
pany, Inc., New York and London, 1945. ix+434 pp. $4.00. 


“This book is designed as a complete elementary exposition of electric-circuit theory requisite in 
the technical foundation of all students of electrical engineering regardless of their expected branch of 
specialization—electric power, communications, or electronics” (from Author's Preface). As such it may 
be recommended to readers of the Quarterly, experts in other than electrical fields, who may at times have 
difficulty in following the exposition of mathematical methods as applied to electrical problems. They 
will find in this book by Professor Frazier, of Massachusetts Institute of Technology, a modern presenta- 
tion of the field of remarkably broad coverage in a relatively small volume. The power and generality of 
modern methods, such for instance as the various types of network transformations, are very well pre- 
sented and thoroughly exemplified. The author has taken great pains to point out possible pitfalls, and 
if his reader will give equally great attention to details he will find himself amply repaid. Historical refer- 


ences and a selected bibliography enhance the value of this book. 
P. LECORBEILLER 


Transmission lines, antennas and wave guides. By Ronold W. P. King, Harry Rowe 
Mimno and Alexander H. Wing. McGraw-Hill Book Company, Inc., New York 


and London, 1945. xv +347 pp. $3.50. 


The book is divided into four chapters. The first chapter, on transmission lines, is written by Alex- 
ander H. Wing; the second and third, respectively on antennas and on wave guides, is by Ronold W. P. 
King; the short concluding chapter is on wave propagation by Harry Rowe Mimno. 

The chapter on transmission lines concentrates on those topics which in recent years have interested 
research workers in microwave laboratories. Those parts of the theory which are needed in problems of 
long line communication, such as crostalk and interference problems, are not considered; but ample 
attention is given to high frequency measurements, impedance matching, suppression of harmonics, etc. 
The emphasis is definitely on high frequencies and on relatively short lines .The exposition is good. 

The chapter on antennas constitutes one-half of the book. For this reason it is particularly unfor- 
tunate that it should contain so much misinformation and misinterpretation. For the most part it would 
be difficult for an inexpert reader to recognize what is right and what is wrong. Throughout, the reader is 
given to understand that the conclusions are based on rigorous electromagnetic theory. Engineering ap- 
proximations in common use are called “very crude” if they are in error by as much as twenty-five per cent 
and one is led to believe that those approximations which are called “good” by the author are really good. 
Apparently, however, the author has not set a uniform objective standard of quality of approximations. 
He declares that his theoretical impedance curves are in “good agreement” with measured impedances. 
He does not give the measured values; but measured values from three published sources, and one un- 
published but made known to the author, agree among themselves and disagree with King’s curves, in 
some regions by as much as twenty-five to seventy per cent. These measured values also agree with the 
theoretical results published by this reviewer and by Marion C. Gray. These facts are not mentioned in 
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the book. The author's attitude seems to be that expressed in one of his latest papers (Jour. App. Phys., 
August, 1945, p. 445): “In many instances disagreement between experimental and theoretical results 
may be a better check on the theory than close agreement.” 

On pp. 104 and 107 are shown curves relating to the length of the antenna at resonance (which is 
defined as the condition for which the input reactance vanishes) and the corresponding input resistance. 
On each curve, there is a point marked “sphere.” The captions explain that the sphere is regarded as a 
cylinder whose height is equal to the diameter. One of these points is taken from a book by J. A. Stratton 
and the other from a paper by E. B. Moullin. The former was calculated for free and not forced oscilla- 
tions; in fact, in the case of a transmitting spherical antenna the input susceptance is always capacitive 
and the input reactance does not vanish. E. B. Moullin calculated an approximate re-radiation resistance 
with reference to the maximum current of a sphere in a certain impressed field and not the input resist- 
ance of the spherical transmitting antenna. In fact, the latter resistance depends markedly on the separa- 
tion between the two halves of the spherical antenna; if this separation is zero as implied by the author of 
the antenna chapter of the book, the input resistance becomes equal to zero automatically. 

At times the author brands a correct conclusion as incorrect and then gives an incorrect result to 
replace it. For example, on p. 223 he purports to show that the effective area of a “half wave” self-tuned 
antenna depends considerably on its radius. He assumes that the effective length of the antenna is inde- 
pendent of the radius and takes into consideration only the variation of the effective area with input 
resistance. Actually, the effective length also varies with the radius and if this effect is included, the effec- 
tive area of the half-wave antenna is found to be nearly independent of the radius—a conclusion well 
known in the art. 

The chapter on wave guides occupies a relatively minor position in the book. It is confined primarily 
to detailed descriptions of various types and modes of propagation and the facts are substantially accu- 
rate. The inequality (10.1) on p. 251 is unduly restricted; but the fault is not particularly serious. On p. 
269 we find: “The upper frequency limit of the 7.Mo,: mode from the point of view of single mode opera- 
tion is the cut-off for the 7E,.,; mode.” The statement is not true; but it is clearly an over-sight and is not 
likely to cause serious trouble. 

The book is concluded with an excellent thumb-nail sketch of factors affecting wave propagation 
over the earth. It is hard, however, to pass without comment the author’s apparent approval of recent 
efforts to ascribe specific meaning to such general terms as “low, medium and high frequencies.” If these 
recommendations are put into effect, the language will needlessly be robbed of valuable general terms. 

S. A. SCHELKUNOFF 


Theory of flight. By Richard von Mises with the collaboration of W. Prager and 
Gustav Kuerti. (McGraw-Hill Publications in Aeronautical Science, Jerome C. 
Hunsaker, Consulting Editor.) McGraw-Hill Book Co., Inc. New York and 
London, 1945. XII+629 pp. $6.00. 


This very comprehensive engineering text book is different from similar books in the same class; the 
author's extensive knowledge of the basic theories and the fundamental principles is everywhere evident. 
According to the preface the book is written primarily for new graduate students. However many of the 
chapters require a most thorough preparation in applied mechanics and considerable insight in fluid 
dynamics. The book will be of considerable interest to engineers who wish to familiarize themselves with 
particular aspects of the problems of engineering aerodynamics. The chapters on airplane performance 
control and stability are particularly complete with numerous useful references to experimental results. 

THEODORE THEODORSEN 


The simple calculation of electrical transients. By G. W. Carter. Cambridge: At the 
University Press, New York: The Macmillan Company, 1945. viii+120 pp. $1.75. 


In this little book Mr. Carter explains how to use Heaviside’s operational method in the transient 
analysis of linear networks consisting of a finite number of meshes and does it very well. The method is 
explained step-by-step and each step is illustrated by practical examples. 

The book is addressed to the engineer who wants to be able to use the operational method with con- 
fidence but is willing to accept some rules on faith. In the brief introductory chapter the reader learns the 
characteristics of the circuits to which he can apply the method. There a parallel is drawn between the 
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Steinmetz method of steady state analysis and the Heaviside method of transient analysis. In Chapter 2 
the differential equation of a simple circuit is solved by successive integrations; this permits the author to 
introduce the operator Q, standing for = which is easier to understand than the Heaviside operator p. 
It is only after some experience with Q that p is brought into the picture. Gradually the method is de- 
veloped and the reader learns to apply it to initially “dead” circuits and then to circuits in any initial 
state. The book begins with very simple examples, and it ends with complicated ones; thus, it should be 


easy for the student to gain confidence in the application of the method. 
S. A. SCHELKUNOFF 
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R.T.P. Translation No. 2468, The mechanics of the plastic deformation of mild steel. 
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Schleicher. 14 pages. 
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